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ABSTRACT 

We have modeled in detail the evolution of rich open star clusters such as the Pleiades, 
Praesepe and Hyades, using simulations that include stellar dynamics as well as the 
effects of stellar evolution. The dynamics is modeled via direct TV-body integration, 
f^S , while the evolution of single stars and binaries is followed through the use of fitting 

formulae and recipes. The feedback of stellar and binary evolution on the dynamical 
| evolution of the stellar system is taken into account self-consistently. 

C*~) • Our model clusters dissolve in the tidal field of the Galaxy in a time span on the 

order of a billion years. The rate of mass loss is rather constant, 1 ~ 2 M per million 
years. The binary fraction at first is nearly constant in time, then increases slowly near 
O ^ the end of a cluster's lifetime. For clusters which are more than about 10 s years old 

the fractions of stars in the form of binaries, giants and collision products in the inner 
few core radii are considerably higher than in the outer regions, beyond the cluster's 
half mass radius. When stars with masses £ 2 Mq escape from the cluster, they tend 
to do so with velocities higher than average. 

The stellar merger rate in our models is roughly one per 30 million years. Most 
mergers are the result of unstable mass transfer in close binaries (~ 70%), but a 
k> \ significant minority are caused by direct encounters between single and binary stars. 

j_j ■ While most collisions occur within the cluster core, even beyond the half mass radius 

C$ ' collisions occasionally take place. We notice a significant birthrate of X-ray binaries, 

most containing a white dwarf as the donor. We also find some X-ray binaries with a 
neutron-star donor, but they are relatively rare. The persistent triple and higher order 
systems formed in our models by dynamical encounters between binaries and single 
stars are not representative for the multiple systems observed in the Galactic disk 
or in open clusters. We conclude that the majority of multiples in the disk probably 
formed when the stars were born, rather than through later dynamical interactions. 

Key words: Methods: N-body simulations - Binaries: general - Binaries: close - 
Open cluster and associations: general - Open cluster and associations: NGC2516, 
NGC2287, Praesepe, Hyades, NGC 2660, NGC 3680 



1 INTRODUCTION 

Open clusters are useful laboratories for studying the inter- 
play between single star evolution, binary star evolution and 
stellar dynamics. Unlike their bigger (and older) siblings, 
the globular clusters, they contain a manageable number of 
stars, and the evolution of the majority of the stars is not 
expected to be strongly affected by the dynamics. Still, a 
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significant number of collisions and subsequent stellar merg- 
ers can take place, as well as dynamically induced exchange 
reactions between single stars and binaries. For all these 
reasons, the usual zoo of objects created in binary stellar 
evolution is significantly enlarged by the presence of even 
more exotic specimens that could not have formed in vitro 
through isolated evolution, but only in vivo through the dy- 
namical interplay of initially unrelated (single or multiple) 
stars. 

The simulations reported here have been run on a 
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Table 1. Observed and derived parameters for several open star clusters with which our simulations may be compared. References to 
the literature (second column) arc: (a) Pinfield ct al. (1998); Raboud & Mermilliod, (1998); Bouvicr ct al (1998); (b) Abt & Levy (1972); 
Dachs, J & Kabus (1989); Hawley et al. (1999). (c) Harris et al. (1993); Ianna ct al (1987); Cox (1954). (d) Andricvsky (1998); Jones & 
Stauffer (1991); Mermilliod & Mayor (1999); Mermilliod ct al. (1990); Hodgkin et al. (1999). (e) Perryman et al. (1998 and references 
therein) Reid & Hawley (1999); (f) Frandsen et al. (1989); Hartwick & Hesscr (1971); Sandrelli et al. (1999). (g) Hawley et al. 1999 ; 
Nordstrom et al. (1997); Nordstrom et al. (1996), Subsequent columns give (3) the distance to the cluster (in pc), (4) the cluster age (in 
Myr), (5) the half mass relaxation time, (6) the crossing time, (7) the total mass (in Mq), (8) estimate for the half mass radius (in pc), 
and (9) the core radius. In cases where the relaxation time is not given in the literature, we calculate it (see Paper IV); these entries are 
indicated by *. The final two columns contain information about cluster stellar content. The column labeled / s: b:t indicates the number 
of single stars, binaries and triples (separated by colons). For clusters where the numbers are given directly by observations, the table 
gives the observed numbers of each system. If the binary fraction is derived by other methods, wc give the relative fractions normalized 
to the number of single stars. The last column (A/b S s:Be:gs) gives the number of observe blue stragglers, Be stars and giants, separated 
by a column. 
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Table 2. Initial conditions and parameters for the selected models. The first column gives the model name, followed by the cluster mass 
(in Mq), the King (1966) parameter Wo, the distance to the Galactic center (in kpc), the initial relaxation time and half mass crossing 
time (both in Myr). The remaining columns give the location of the cluster's Jacobi surface along the X— (towards the Galactic center), 
Y — and Z— (towards the Galactic pole) axes, and the initial virial radius, half mass radius and core radius (all in parsec). 
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GRAPE-4 (Makino et al. 1997) system, a special-purpose 
computer designed to speed up stellar-dynamical calcula- 
tions. While the models in this paper contain only ~ 3,000 
stars, appropriate for open clusters, we have started to use 
the next-generation special-purpose computer, the GRAPE- 
6, to extend our simulations to include one or two orders of 
magnitude more stars. This will enable us to model glob- 
ular clusters on a star-by-star basis, including those with 
very dense cores where most stars are influenced by various 
'traffic accidents' in the form of encounters and mergers. 
The simulations reported here thus play a double role: as- 
trophysically, they provide insight in the evolution of open 
clusters, and computationally, they help pave the way for 
the modeling of the more complex globular clusters. 

In this paper we describe a series of self-consistent sim- 
ulations of the evolution of open star clusters. All important 
effects are included to some degree of realism: stellar evolu- 
tion, binary evolution, the internal dynamical evolution and 
the effects of the tidal field of the Galaxy. This is the fifth 
paper in a series in which we have gradually increased the 
'ecological' complexity of stellar interactions to a realistic 
level. In this paper we analyze the results of the same eight 
calculations performed for our previous Paper IV (Porte- 
gies Zwart et al. 1999). There we concentrated on global 



cluster properties and compared them with observed open 
clusters, such as the Pleiades, Praesepe and Hyades. Specifi- 
cally, we took a photometric standpoint and studied changes 
in the Hertzsprung- Russell diagram and cluster morphology 
as functions of time and initial conditions. Here we concen- 
trate on the binary population and on higher-order multiple 
systems, and we study the dynamical and observational ef- 
fects of these binaries and multiples on the evolution of the 
stellar system as a whole. 

Detailed descriptions of the numerical methods used 
and global assumptions made in our calculations are given 
in Paper IV. Since we analyze the same data from a different 
perspective we do not repeat this information here. Instead, 
we quickly review the initial conditions of our models, and 
restate the methods used; for more details, see Paper IV. 



2 METHODS 

In this section we briefly discuss the initial conditions of 
our models and the numerical techniques used in our model 
calculations. For details, refer to Appendix B and C of Paper 
IV (Portegies Zwart et al. 1999). 
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Table 3. Initial conditions for the stellar and binary population. 
The first column gives the parameter, the second column gives 
the functional dependence, followed by the lower and upper limits 
adopted. (RLOF indicates that the initial binary may be semi- 
detached.) 
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solar neighborhood, we find that a model with Wo = 6 has 
to be placed slightly farther (10.4 kpc) from the Galactic 
center than the Sun, while a model with Wo = 4 has to 
reside somewhat closer (6.8 kpc), in order to obey the above 
relationships. 

For a total cluster mass of 1600 Mq, the Lagrange points 
of our two standard clusters lie at distances 14.5 pc (Wo = 
4) and 21.6 pc (Wo = 6), from the cluster center. Stars 
are removed from a simulation when their distance to the 
cluster's density center exceeds twice the distance to the first 
Lagrangian point. Tabled reviews the selected parameters 
and initial conditions. 



2.1 Initial conditions 

We are interested in moderately rich (~ 1000 stars) open 
clusters of intermediate age ( ,$ lGyr). Starting from the 
currently observed mass and dynamical properties of such 
clusters we have reconstructed two plausible sets of initial 
conditions, which are presented in Table|!2] Each calculation 
is performed four times with different random seeds in order 
to improve statistics. The notation in this paper is identi- 
cal to that used in Paper IV (see its Appendix A for an 
overview) . 

For our simulations we assume a Scalo (1986) initial 
mass function, with minimum and maximum masses of 
0.1 M and 100 Mq, respectively, and mean mass (m) ~ 
O.6M0 at birth. Consistent with the above mass estimates, 
our simulations are performed with 1024 single stars and 
1024 binaries, for a total of 3096 (3k) stars. 

Stars and binaries within our model are initialized as 
follows. A total of 2k single stars are selected from the initial 
mass function and placed in an equilibrium configuration in 
the selected density distribution (see below). We then ran- 
domly select lk stars and add a second companion star to 
them. The masses of the companions are selected randomly 
from the IMF between 0.1 M and the primary mass (see 
Tab.|3J. For the adopted binary fraction, the restricted sec- 
ondary mass range translates into an overall mean mass of 
O.53M . Then the other binary parameters are determined. 
Binary eccentricities are selected from a thermal distribution 
between and 1. Orbital separations a are selected with uni- 
form probability in log a between Roche-lobe contact and an 
upper limit of 5,000 A.U. (= 10 6 R Q , or 0.02pc). When a 
binary appears to be in contact at pericenter, that particu- 
lar orbit choice is rejected and new orbital parameters are 
selected. Table|3] gives an overview of the various distribu- 
tion functions from which stars and binaries are initialized. 
(Figure [I] shows the initial distributions of orbital period 
and eccentricity.) 

We select initial density profiles from the density dis- 
tributions described by Heggie & Ramamani (1995) with 
Wo = 4 and Wo = 6, and refer to these models as W4 
and W6, respectively, throughout this paper. We have per- 
formed four independent simulations for each set of initial 
conditions, labeled / to IV. 

All our cluster models start with the same virial radius 
of Ro — 2.5 pc. This implies a conveniently constant scaling 

between the cluster dynamical time scale ~ (GMo/i?o) ^ 2 
(= 1.5 Myr for M = 16OOM ) and the time scale for stellar 
evolution. Each cluster is assumed to precisely fill its Jacobi 
surface at birth. Given the observed Oort constants in the 



2.2 The JV-body integrator 

The TV-body integration in our simulations is carried out 
using the kira integrator, operating within the Starlab 
software environment (McMillan & Hut 1996; Portegies 
Zwart et al. 1999) . Kira uses a fourth-order Hermite scheme 
(Makino & Aarseth 1992), incorporates block time steps 
(McMillan 1986a; 1986b; Makino 1991), includes special 
treatments of close two-body and multiple encounters of ar- 
bitrary complexity, and a robust treatment of stellar and 
binary evolution and stellar collisions (see below). As men- 
tioned above, the special-purpose GRAPE-4 (Makino et al. 
1997) system is used to accelerate the computation of grav- 
itational forces between stars. 

2.3 Stellar evolution 

The evolution of the stars is adapted from the prescription 
by Portegies Zwart & Verbunt (1996, §2.1). However, some 
changes are made to the mass loss in the main-sequence 
stage for massive stars and for mass-transfer remnants. We 
follow the details of model B from Portegies Zwart & Yun- 
gelson (1998). For our treatment of stellar mass loss, see 
Portegies Zwart et al. (1998). 

Neutron stars receive a high-velocity kick at the mo- 
ment of their formation. The reasons for the occurrence of 
these high kick velocities have been discussed in detail by 
Portegies Zwart & van den Heuvel (1999) and the implica- 
tions for binary evolution are discussed by Portegies Zwart 
(2000). The distribution from which the kick velocity should 
be taken is less certain. We adopt the velocity distribution 
proposed by Harmann (1997) which has a dispersion veloc- 
ity of 450 km s . Each neutron star formed received a kick 
chosen from this distribution and oriented in a random di- 
rection. 



3 RESULTS 

3.1 Overall evolution of the models 

In this section we review the global properties of models 
W4 and W6, summarizing the more extensive discussion 
presented in Paper IV. Where we discuss individual mod- 
els, we focus on the same two models discussed in Paper 
IV (model W4-IV and W6-III). In other cases the data for 
several models, or even all models, are combined in order 
to improve statistics; those cases are indicated specifically 
below. 
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Figure 1. Total mass as a function of time for model W4 (dashes) 
and W6 (solid). 

Figures and [5] present the overall evolution for our 
two sets of cluster initial conditions. Figure Q shows the 
time evolution of the cluster mass. Models W4 and W6 lose 
mass at roughly constant rates of about 1.4 M© (0.09%) and 
0.82 Mq (0.05%) per million years, respectively. These mass 
loss rates result in the disruption of the cluster at about 1 
Gyr for model W4 and around 2 Gyr for the more concen- 
trated model W6. The higher mass loss rate of model W4 
is mainly attributable to its closer proximity to the Galac- 
tic center and not per se to its lower concentration. These 
rates are consistent with a mass loss rate of 1.2 Mq per mil- 
lion years for the slightly more massive star cluster studied 
by Hurley et al. (2001; 2002). The mass-loss rates per half 
mass relaxation time derived by Portegies Zwart et al. (2001) 
for dense star clusters near the Galactic center are about a 
factor of four higher than the corresponding rates for the 
clusters studied here. 

Figure |5] shows the evolution of the Lagrangian radii 
containing 5%, 25%, 50% and 75% of the cluster mass, 
for models W4 and W6. Both clusters expand by about 
a factor of three during the first half-mass relaxation time 
(~ 140 Myr). It is not clear whether the clusters experience 
any significant core collapse during this early phase. In Pa- 
per IV we concluded that the clusters do not experience core 
collapse due to the effects of binary activity, which tends to 
counteract core contraction (see McMillan et al. 1990, 1991). 
At later times (t ^ 650 Myr for model W4 and somewhat 
later for model W6), the Lagrangian radii decrease again as 
the cluster starts to dissolves in the tidal field of the Galaxy. 



3.2 Global binary properties 

Both models start with 50% binaries, so two-thirds of the 
stars are initially members of binary systems. The adopted 
upper limit (10 6 Kq) on the initial semi-major axis means 
that most initial binaries are relatively wide. The fraction of 
hard (\E\ > IkT) binaries was 0.46 ±0.01 for all models, ir- 
respective of the initial density profile. (The thermodynamic 
unit kT is defined in Paper IV). 

Figure [3] shows the evolution of the binary fraction. 
Shortly after the start of the runs, the binary fractions drop 
to about 42% and 48% for models W4 and W6, respectively, 
and remain roughly constant thereafter. This initial decrease 
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Figure 2. Evolution of the 5%, 25%, 50% and 75% Lagrangian 
radii for model W6 (solid lines for the 5% and 50%, dashed lines 
for the 25% and 75% Lagrangian radii). Lagrangian radii for 
model W4 are presented as dotted lines. 
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Figure 3. Binary fractions as a function of time for models W4 
(dotted line) and W6 (solid line). 



in binary frequency is consistent with the disruption of all 
soft (\E\ 5, IkT) binaries. Supernova explosions dissociated 
1.8% of the binaries in the W4 models, and about half that 
number in the W6 models. This difference is mainly the 
result of random fluctuations. Fluctuations in the binary 
fraction after ~ 100 Myr are mainly the result of escapers 
and mergers. The number of mergers resulting from unstable 
mass transfer or direct collisions is 12.1 ± 0.5% for models 
W4 and W6. The escape rate of binaries closely follows the 
escape rate of single stars; both are relatively constant with 
time. Model W6 loses one binary per 2Myr; the escape rate 
for model W4 is about twice as high. 



3.3 Evolution of binary parameters 

The ecological interplay between stellar (and binary) evolu- 
tion, stellar (and binary) dynamics, and the external tidal 
influence of the Galaxy transforms the initial distributions 
of stars and binaries in complex ways. We assume that ini- 
tially all binaries are distributed throughout the cluster in 
the same way as single stars. However, since binaries are on 
average more massive than single stars, they subsequently 
tend to sink to the cluster center, where superelastic encoun- 
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ters quickly modify the spatial distributions of both stars 
and binaries. 

These encounters come in many forms, from the more 
common single-star-binary and binary-binary types to the 
rare forms that involve triples and occasionally even higher- 
order multiple systems. The combined effects of such en- 
counters significantly modifies the binary distribution func- 
tions, in terms of radius and binding energy as well as stel- 
lar type. In addition, stellar evolution and isolated binary 
evolution also cause binary parameters to evolve in time, in- 
dependent of the occurrence of dynamical close encounters. 
In this subsection, we will present some of the net results of 
these complex processes. 



3.3.1 Distributions of orbital period and eccentricity 

Figure shows the initial distribution of orbital period (in 
years) and eccentricity of a subset of the primordial bina- 
ries of model W6 (the initial conditions for model W4 are 
identical). Figures 2t> and c show the same distributions for 
models W4 (panel b) and W6 (panel c) at 600 Myr. In the 
middle panel, for model W4, only 871 binaries remained. In 
order to facilitate visual comparison, the upper and lower 
panels also show only 871 binaries out of the larger numbers 
that could have been plotted. 

A striking feature of the three panels in Fig|l] is the lack 
of wide (p or b ^ 10 4 yr) binaries at later times. These binaries 
are completely absent in model W4, while in model W6 a 
small population of wide binaries remains. This deficiency 
in wide binaries in evolved clusters is caused by ionization of 
the soft binaries by dynamical encounters (see also Figure 
0J. For model W4 this process is more efficient than for 
the more concentrated models because 1) the W4 models 
experience a high density phase during their early evolution, 
and 2) the W6 models have extended outer regions with 
relatively low stellar densities because of the weaker tidal 
field. All surviving binaries with p or b Si 10 year are located 
well beyond the cluster's half mass radius, and most contain 
at least one white dwarf. In many cases, the orbital periods 
of these non-interacting binaries has increased substantially 
due to the large amount of mass lost from one (or both) of 
the component stars. 

Binaries located above and to the left of the left-most 
dotted curves experience strong tidal effects during their 
early evolution. As a result, they quickly circularize, as can 
be seen by their precipitation to the zero-eccentricity line 
at the bottom in the last two panels. The majority of these 
systems contain at least one white dwarf or helium star. Oc- 
casionally a binary may enter the 'forbidden' left-most area 
at a later time as a result of a strong dynamical encounter. 
Such a binary will either be quickly circularized or its com- 
ponents will collide and merge to form a single star (see 

MM - 

Figure shows how the percentage of tidally circular- 
ized binaries evolves over time. Not surprisingly, the very 
hard binaries, which have high binding energy and therefore 
tight orbits, have on average undergone much stronger tidal 
interactions and therefore have had a much larger chance to 
circularize. The four percentages presented in the figure — 
hard binaries and all binaries for each of the two classes 
of model — at first change little during the evolution of the 
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Figure 4. Scatter plots for the orbital period versus the eccen- 
tricity for 871 binaries for various models and moments in time: 
a) initial distribution, b) the distribution for model W4 at an 
age of 600 Myr, c) model W6 at t = 600 Myr. Each panel is con- 
structed using 871 binaries. Contours give the 0.5, 0.25 and 0.125 
iso-density curves of the distribution in the two coordinates plot- 
ted here. The left dotted curve indicates at which orbital period 
and eccentricity a binary with a 0.1 Mq zero-age primary circu- 
larizes during the time span of our simulations. The right dotted 
curve is for a 2.26 Mq primary at zero age. 

cluster, but they increase at later times in the case of model 
W4. 

Note that in Fig^] there are some circularized binaries 
with long orbital periods, far to the left of the right-most 
dotted line, in the middle and lower panel. These are absent 
in the upper panel, which shows the initial conditions. The 
tidal circularization of these systems occurred when one of 
the stars ascended the giant branch before the onset of mass 
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Figure 5. The relative number of binaries that are tidally cir- 
cularized, as a function of time, for the models W4 (dashed lines) 
and W6 (solid lines). In each case, the lower lines present the over- 
all fraction of binaries that are circularized, while the top curves 
indicate the fraction of circularized binaries among those binaries 
that are very hard (\E hin \ > 1000/cT). 




ecc 

Figure 6. Cumulative distribution of the binary eccentricity. 
Solid line, initial conditions (the thermal distribution). The 
dashed and dotted line represent the eccentricity distributions 
at an age of 600 Myr of binaries with a binding energy Bbin < 
-10 3 A:T and 

^bin < — 10 4 fcT, respectively. For the latter two 
curves the data for all models W4 and W6 are combined. How- 
ever, we excluded those binaries which have been circularized by 
tidal forces. 



transfer. While in many cases this type of evolution leads to 
a later shrinking of the orbit, in some rare cases the binaries 
remain wide after a period of modest mass transfer. 

Since tidal forces drop off quickly with distance, in most 
cases they either succeed in fully circularizing a binary, or 
they do not cause much of an effect. Consequently, the eccen- 
tricity distribution of most non-circularized binaries changes 
little with time. Figure |5] illustrates this by presenting the 
cumulative distribution of all binaries with the exception of 
those with zero eccentricity. Note that this general argument 
does not hold for the hardest binaries: when we select only 
the hard (\E hin \ > 1000/fcT) or super-hard (|-E bIn | > 10 4 fcT) 
binaries, clear deviations from the initial thermal distribu- 
tion are evident. This is not surprising: a 10 4 fcT binary sim- 
ply does not have enough room to develop an eccentricity 
of, say, 0.9, because the semimajor axis is too small. 



3.3.2 Radial distribution of binaries 

FigureQshows the binary fraction as function of distance to 
the cluster center at birth and at t = 600 Myr, adopted as 
a representative snapshot of an older cluster. (Qualitatively 
similar results would have been obtained had we selected 
any other time between ~ 200 Myr and 1 Gyr.) The upper 
solid line shows the initial distribution for all binaries, re- 
gardless of their binding energy. The lower solid line presents 
the initial fraction of binaries with a binding energy of at 
least |i?bin| > lOOOfcT. To avoid clutter, we have shown the 
statistical uncertainties in the form of error bars only for the 
bullet points; the uncertainties in the open and close trian- 
gles can be estimated by the scatter between neighboring 
points. 

The triangles in Figure |7| shows the distribution for all 
binaries at t — 600 Myr for model W6 (open) and for model 
W4 (filled). The two distributions show similar trends, ex- 
cept for a generally smaller binary fraction in model W4 (see 
also Figure |5J. The binary fraction in the central portion of 
the cluster is considerably higher than it was initially (upper 
solid line), because of mass segregation. For higher r values, 
the fraction of objects that are binary remains below the 
original value, all the way to and beyond the tidal radius. 

The bullets indicate the distribution of binaries with 
l-Ebinl > 1000/cT at a cluster age of 600 Myr. In other to 
improve our statistics, we have combined the data from the 
two models; the distributions are similar, with model W6 
containing ~ 8 % more hard binaries than model W4. Note 
that we find a higher proportion of hard binaries over the 
entire cluster field as compared to the initial distribution. 
This is mainly a result of dynamical evolution: hard bina- 
ries tend to get harder, and therefore some of the binaries 
contributing to the bullet points started off with a binding 
energy less that lOOOfcT, and thus were not counted in the 
construction of the lower solid line. 

The higher central concentration of very hard binaries, 
as compared to less hard binaries, is evident in Figure |HJ 
which figure shows the cumulative distribution for single 
stars, hard |_B b in| > 1000/fcT and super hard \E hin \ > W 5 kT 
binaries. The figure also presents the cumulative radial dis- 
tribution for higher order (mostly triple) systems. The left- 
most and right-most dotted curves in Figure |S] show the 
radial distribution for model W6 at birth and at an age of 
600 Myr, respectively. Initially, both single stars and bina- 
ries follow the left-most dotted curve. At an age of 600 Myr, 
all binaries are more centrally concentrated than the single 
stars (see the figure caption). The harder the binaries, the 
more centrally concentrated their distribution has become. 
In the case of triples, the condensation toward the center is 
extreme, with a triple half-mass radius of around 0.3 pc. 

The high central concentration of the very hard |£bm| > 
10 5 fcT binaries is in part due to their early history. These 
binaries are largely the product of an unstable phase of 
mass transfer. The stellar components in these binaries were 
therefore more massive than they are at t = 600 Myr, so they 
are more more strongly affected by mass segregation, causing 
them to sink to the central regions. In contrast, the mod- 
erately hard binaries with binding energies between lOOfcT 
and 10kT (not shown in Figure |SJ are barely more centrally 
concentrated than the single stars. The triples are strongly 
centrally concentrated. This is a the most dramatic result 
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Figure 7. Fraction of binaries N^i n / (N SB + iVbin) as function 
of the distance to the cluster center. The upper solid horizontal 
line shows the initial binary fraction over the entire star cluster. 
The cut off at about 10 pc is close to the location of the tidal 
radius. (Initially there were no stars beyond the tidal radius.) 
The lower solid curve presents the initial distribution for binaries 
with |i?binl > lOOOfcT. The open triangles give the distribution of 
all binaries for model W6 at an age of 600 Myr. The solid triangles 
present the distribution of all binaries for model W4 at an age 
of 600 Myr. The bullets (with 1 a Poissonian error bars) present 
the fraction of hard |i?bin| > lOOOfcT binaries as a function of the 
distance to the cluster center at an age of 600 Myr. These data 
combine models W4 and W6. To guide the eye we have placed a 
straight (dotted) line through these bullet points. 




r [pc] 

Figure 8. Cumulative radial distribution of single stars and bina- 
ries in models W4 and W6 (combined). The cumulative distribu- 
tion for all single stars and binaries at birth and at t = 600 Myr 
are given by the left and right dotted curves, respectively. The 
solid curve gives the distribution for binaries with |i?bml > 10 3 fcT 
at an age of 600 Myr. The dashed curve gives the same distribu- 
tion for |-Ebinl > 10 5 fcT binaries. The dash-3-dotted curve gives 
the radial distribution for higher order (mostly triple) systems. 
To improve statistics we accumulated all triples for models W4 
and W6 over the time range of 550 Myr to 650 Myr. 



of dynamical interactions: triple systems were not initially 
present, and must be created dynamically through encoun- 
ters. Since the encounter rate increases sharply toward the 
cluster center, triples and higher-order systems are born, 
and often quickly destroyed again, near the very center of 
the cluster. 
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Figure 9. Cumulative distribution of the velocities of escaping 
single stars and binaries in models W4 and model W6 (combined). 
The two dotted lines show the velocity distributions of single stars 
in the models at zero age (right) and at t = 600 Myr (left). The left 
solid curve presents the distribution of the escaper velocity of all 
stars and binaries integrated over time. The right solid curve gives 
the velocity distribution of stars which escaped within the first 
half-mass relaxation time (~ 150 Myr). The two dashed curves 
give the escape speeds for stars with masses m CBC > 2 Mq (upper) 
and m CBC > 4Mq (lower). 



3.4 Escapers 

The most important mechanism by which stars escape from 
the cluster is tidal stripping. Stars are assumed to have been 
tidally stripped once they reach two Jacobi (tidal) radii from 
the cluster center; they are then removed from the simula- 
tion. Stripped stars generally leave the cluster with rela- 
tively low velocities, as illustrated in Figure [5] which shows 
the distribution of escaper speeds for models W4 and W6. 

It is convenient to draw a distinction between tidal 
evaporation and dynamical ejection of stars. Operationally, 
we distinguish between the two processes by the speed of the 
escaper; an "ejected" single star or binary leaves the cluster 
with a speed exceeding the escape velocity. That is, we in- 
fer the dynamics from the speed. Typically, this velocity is 
imparted to the escaping object during a strong interaction 
with another object in the cluster core. Because of the rela- 
tively low central densities (and hence low reaction rates) in 
our model clusters, this type of ejection is rather uncommon. 
White dwarf ejection velocities are not significantly differ- 
ent from those of main sequence stars, but giants have on 
average somewhat lower escaper velocities. This is mainly a 
consequence of the fact that giants are generally large and 
therefore cannot have very close encounters without collid- 
ing, so they are less able to pick up large escape velocities. 

Another important escape mechanism is supernova ex- 
plosions. A neutron star formed in a core-collapse supernova 
is hurled into space with a "kick" velocity that can easily ex- 
ceed several hundreds of kilometers per second, greater than 
the cluster escape speed by 1-2 order of magnitude. These 
objects are hidden in the high-velocity tail of Figure|5] how- 
ever, they are clearly visible in Figure [THl as the population 
of high-speed single objects with masses between 1 Mq and 
2M Q . 

A supernova in a binary system has several effects on 
the velocity distribution of both single and binary stars. 
The sudden mass loss in the supernova event, as well as 
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the asymmetric velocity kick, affect a binary's internal or- 
bital parameters and its center of mass velocity. When the 
binary is dissociated in the supernova, the companion to the 
exploding star is ejected with its orbital velocity, while the 
newly formed compact object receives an additional velocity 
kick. The effects of supernovas on binary systems, and on 
the velocities of the binaries and single stars which result, 
are reviewed by Portegies Zwart (2000). 

Figure [§] shows the velocity distribution of escaping 
stars and binaries from models W4 and W6. Due to the 
similarities in escape speeds between the two models, we 
have combined the data in a single plot. The escaper veloc- 
ity distribution for single stars is slightly higher than that 
of the binaries. For clarity we show only the combined dis- 
tribution of single and binary escapers in Figure [5] Escaper 
velocities are generally only slightly higher than the cluster 
velocity dispersion, and somewhat higher in the W4 models 
than in W6. The small fraction of high-velocity escapers in 
Figure [5] is due to neutron stars escaping after supernovae 
explosions. 

In addition to the overall escaper velocities for the com- 
bined models, Figure|J|]also shows the distribution for stars 
which escaped within the first half-mass relaxation time 
(right solid line). This distribution has a considerably higher 
mean than the overall distribution. The figure also shows 
the velocity distributions for escaping stars more massive 
than 2 Mq (top dashed curve) and m csc > 4 Mq (lower 
dashed). These distributions are very different from the av- 
erage escaper speed. Both models W4 and W6 experience 
high-density phases during their early evolution, within 1 
initial half-mass relaxation time. The early escapers there- 
fore have higher velocities. A similar process operates for the 
more massive escapers. The turn-off mass of a 2 Mq star is 
about 800 Myr, while a 4Mq star lives for less than about 
140 Myr. The distribution of escaper speeds for stars which 
escaped within t r i x includes the stars with m > 4Mq. Still, 
the distributions are considerably different, in the sense that 
the massive stars have relatively high escape speeds. The 
reason for this discrepancy is the greater dynamical activity 
of the high-mass stars, which take part much more often in 
dynamical encounters with binaries. 

The dependence of escaper velocity on escaper mass 
is illustrated in Figure 1101 Most escapers have low masses 
and low velocities. The average escaper velocity in the W4 
and W6 runs was 1.14 km s -1 and 1.00 kms, respectively. 
Higher mass stars (and binaries) have considerably higher 
velocities. This trend was first observed by Blaauw (1961, 
see also Gies & Bolton 1986). Binaries are mostly ejected by 
dynamical encounters. Supernovae are also responsible for 
numerous high velocity escapers, both from kicks and from 
binary effects, as mentioned above. The neutron stars are 
clearly visible in Figure [TU1 as the objects having masses of 
about 1.4 Mq and velocities up to ~ 1000 km s" 1 . 

Two rather low-mass (m ~ 0.11 Mq and m~ 0.7 Mq) 
stars have unusually large (~ 100km s _1 and 185km s _1 , 
respectively) escape velocities. The more massive of the two 
was ejected from a tight binary at about t = 5.6 Myr when 
its companion exploded in a Type lb supernova, dissociat- 
ing the binary. The low-mass object was ejected following a 
three-body encounter at t = 223 Myr. 

Of all stars ejected from models W4 and W6, only 109 
(~0.93% out of a total of roughly 12,000) have velocities 
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Figure 10. Mass versus velocity for escaping single stars (dots) 
and binaries (plus signs) in model W6. For the binaries, the total 
mass is used. 



exceeding three times the dispersion velocity i>di sp of the 
parent cluster. Of these, 46 are neutron stars. However, 57% 
of escapers with masses between 4Mq and 11 Mq, and 75% 
of escapers more massive than 11 Mq, have velocities more 
than 3udis P - The average escaper velocity of stars between 
4Mq and 11 Mq is 8.5km s -1 ; stars with masses exceeding 
11 Mq have an average escaper velocity of 21 km s _1 , which 
is much higher than the cluster dispersion velocity. 

These numbers are unusually high compared to the 
number of runaways in the Galactic disk. However, as our 
criterion here we have adopted three times the cluster ve- 
locity dispersion, instead of Blaauw's (1961) criterion of 
40 km s -1 (or three times the velocity dispersion in the 
Galactic disk). The fraction of high velocity stars following 
Blaauw (1961) was between ~ 1 % for stars with m 11 Mq, 
about 2.5% for early B stars and about 20% among O stars 
(m ^ 16 Mq, see also Sone 1991). Among stars of spectral 
type A, Stone (1991) found a runaway frequency of ,$ 0.3 %. 
These numbers are consistent with binary population syn- 
thesis studies of binaries in which stellar dynamical encoun- 
ters are not taken into account (Portegies Zwart 2001). 

The velocity dispersion of all stars in the cluster at 
an age of 600 Myr is «dis P = 0.72 ± 0.23 km s _1 . For stars 
within 1 pc (about the core radius) of the cluster center, 
^disp = 0.98 ± 0.34 km s _1 ; stars between 10 and 15 pc of 
the center have v^ sp = 0.59 ± 0.16 km s^" 1 . Most striking is 
the higher velocity dispersion for stars well outside the tidal 
radius (r > 15 pc) which have «dis P = 0.75 ± 0.20 km s _1 . A 
KS test indicates that the various velocity distributions are 
significantly different. 

Drukier et al. (1998) reported a similar effect in the 
Globular cluster M15, and attributed it to tidal shocking. 
However, this process is not included in our models; we 
therefore conclude that the observed increase in the veloc- 
ity dispersion outside the tidal radius cannot simply be a 
consequence of tidal shocks. 



3.5 Collisions and Coalescence 

Stellar mergers result from unstable mass transfer in a bi- 
nary system or from direct collisions between stars. Some- 
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Table 4. Mergers occurring in models W4 and W6, compared 
with those in model calculations published previously. The first 
column identifies the two stars involved in the merger. The next 
two columns give the number of mergers in each model calcula- 
tion, followed by the fraction of each merger type (expressed as 
a percentage, combining all runs). The last two columns give the 
fractions of mergers which occurred in models S of Papers I and 
II. Mergers between two remnants in Paper II (model S) are in- 
cluded here under {wd, wd}, although a few actually involved a 
neutron star or black hole. 





W4 


W6 


total 


Paper IS 


Paper IIS 


Ntot: 


117 


121 




■ [%} 




{ms, ms} 


92 


110 


84 


77 


80 


{ms, gs} 


19 


11 


12 


14 


6 


{ms, wd} 





1 





5 


11 


{ms, ns} 


1 








1 


1 


{gs, ns} 











3 





{wd, wd} 


5 





2 





2 


{wd, ns} 





2 


1 









times (unstable) mass transfer is initiated by the presence 
of a third star. 

Table^Jgives an overview of the 241 collisions occurring 
during the model calculations discussed here. A schematic 
diagram of the relative frequencies of various collision con- 
figurations in presented in Figure 1111 This figure may be 
compared directly with Figure 2 of Paper I, and with Figure 
5 of Paper II. For easier comparison, we have added to Ta- 
bic 0] two columns summarizing the results of Papers I and 
II. Note, however, that the comparison with papers I and II 
is not really appropriate, as these model calculations lasted 
for 10 Gyr. According to a KS test there is no significant 
difference between the time histories of the collisions in the 
W4 and the W6 models. 

Most ( ^ 80%) mergers in models W4 and W6 occur 
between two main sequence stars (see Table This was 
also the case in Papers I and II. In Paper I, however, only 
collisions between single stars were considered, and the high 
proportion of main-sequence collisions simply reflects the 
high proportion of main-sequence stars. Paper II included 
both collisions between single stars and mergers resulting 
from binary mass transfer, so the results in that case should 
compare better with models W4 and W6. The collisions re- 
ported in this paper, however, have a rather small contri- 
bution from {ms, wd} mergers compared to those in Pa- 
pers I and II. This may be due in part to the smaller ages 
( ^ 1 Gyr) of our models compared with Papers I and II 
(<10Gyr). 

Although clearly limited by small number statistics, an 
apparent difference between models W4 and W6 is the num- 
ber of {wd, wd} mergers, which are much more common in 
model W4. This difference is reflected in the higher Type la 
supernova rate in model W4. Interestingly, the calculations 
in Paper II also had a high {wd, wd} merger fraction. In 
fact, if the binary population of models W4 and W6 were 
allowed to evolve after the disruption of the cluster, the rate 
of {wd, wd} mergers would be very similar to the results of 
Paper II. The cause of the high proportion of white-dwarf 
mergers in model W4 and Paper II is the larger amount of 
dynamical activity in these models, which drives more bi- 
naries into a state of mass transfer, favoring the formation 
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Figure 12. Cumulative distribution of the distance from the clus- 
ter center at which collisions occur in models W4 (dashes) and 
W6 (solid). The dotted line gives the initial distribution of stars 
in the W6 models. 

of white-dwarf binaries (see also Hurley & Shara 2002). The 
resulting Type la supernova rate is discussed in the next 
section. 

Figure [T2l shows the radial distribution of collisions in 
models W4 and W6. More than 80% of all collisions occur 
within the initial half-mass radius (~ 3pc). Collisions in the 
W6 models tend to be more concentrated to the cluster core 
than in the W4 models. It is striking that some (~ 5%) 
collisions occur far ( ^ 8 pc) from the cluster center. These 
collisions are induced by binary evolution, whereas mergers 
in the core are mainly caused by stellar encounters. 

The expected distribution of masses and radii of col- 
liding stars can be calculated following Portegies Zwart et 
al. (1997; see their Eq. 14). Assuming a Maxwellian velocity 
distribution with velocity dispersion v and including the ef- 
fects of gravitational focusing, the number of collisions in 
the cluster per 10 8 year may be expressed as 

where m is the stellar mass. This rate is averaged over all 
other masses. (Note that, based on this estimate, we expect 
no collisions during the time span of our simulations for our 
adopted cluster parameters, underscoring the importance of 
binary interactions.) 

Figure Hoi shows the expected (upper panel; see paper 
III for details) and observed (lower panel) distributions of 
primary and secondary masses of stars involved in collisions 
in models W4 and W6. We excluded here the stars which 
merged as a result of unstable mass transfer. The upper 
panel is created via Eq.^from the initial mass function and 
the same mass-radius relation for zero-age main-sequence 
stars as was used in the model calculations. Gray shades 
indicate collision probability; darker shades correspond to 
higher values. The collisions observed in our models are pre- 
sented in the lower panel. The masses of the colliding com- 
ponents in our model calculations are on average somewhat 
higher than expected. This phenomenon was first observed 
by Portegies Zwart et al. (1999; see also Portegies Zwart & 
McMillan 2002) in simulations of young star clusters having 
high central stellar densities. Portegies Zwart et al (1999) 
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Figure 11. Schematic overview of the relative frequencies of mergers between various stellar types in models W4 and W6 (combined). 
Stellar types are denoted by ms for main sequence, gs for (sub)giants and rm for white dwarfs and neutron stars. Table HI gives the 
numbers of occurrences. 
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Figure 13. Relative collision rates between a primary and a sec- 
ondary star as expected from Eq.0 (upper panel) , and as found 
in models W4 and W6 (lower panel). For the lower panel the 
mergers which resulted from an unstable phase of mass transfer 
are excluded; only the true collisions are taken into account. The 
shading is linear in the encounter probability in the upper panel 
and in the number of collisions in the lower panel. Darker shades 
indicate higher collision frequency. Since the probability distribu- 
tion is symmetric about the line of equal masses, only the lower 
half of each figure is displayed. 



used R 136, a compact and very dense star cluster in the 
30 Doradus region, as template for their simulations. We re- 
turn to this comparison in the discussion section. 

For 240 of the 241 stellar pairs experiencing a collision, 
we were able to reconstruct the orbit before the collision oc- 
curred. Figure [T4l shows the orbital parameters for the bound 
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Figure 14. Orbital periods and eccentricities of coalescing bina- 
ries in models W4 and W6. Notice the rather large population of 
circular (zero eccentricity) binaries which experience a collision. 



systems thus obtained. The single pair for which this recon- 
struction was not possible experienced a collision during a 
supernova event — the neutron star formed in the supernova 
was ejected directly into its 1.9 Mq main-sequence compan- 
ion. In 181 of the 240 collisions (~ 71%), the orbit was cir- 
cular on merger. The eccentricity distribution in the other 
bound cases is consistent with a thermal distribution. This 
result is consistent with the findings in Paper II, where for 
model S, 19% of the collisions were the result of a dynamical 
encounter. 



3.6 Supernovae 

Supernovae can dramatically affect the evolution of a star 
cluster. During a supernova, the exploding star ejects a large 
fraction of its mass at high speed. This mass quickly leaves 
the cluster, as the ejection velocity far exceeds the cluster's 
escape speed. In addition, the newborn compact object may 
also receive a kick velocity large enough to escape the clus- 
ter. Binaries are very fragile to supernovae, and each binary 
hosting a supernova is likely to be dissociated. 



3. 6. 1 Supernova types 

Tables|^|and|S|present overviews of the supernovae observed 
in the models discussed in this paper, and compare them 
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Table 5. Observed supernova types (Obs, from Cappellaro et al. 
1997) and supernovae occurring in a standard Scalo (1986) mass 
function (model IMF), the population synthesis calculations of 
model AK by Portegies Zwart & Verbunt (1996, model PZV-AK) 
and in models W4 and W6, and when the dynamical evolution 
of the stellar system is ignored (models W4-nd and W6-nd). The 
first column indicates the model, followed by the number of su- 
pernovae of types la, lb, Ic and II. All numbers are normalized to 
1. 



model 


la 


lb 


Ic 


II 


Obs 


0.09 


0.11 


0.18 


0.62 


IMF 






0.12 


0.88 


PZV-AK 


0.03 


0.13 


0.12 


0.72 


W4 


0.22 


0.20 


0.07 


0.50 


W4-nd 


0.06 


0.20 


0.17 


0.57 


W6 


0.10 


0.20 


0.10 


0.63 


W6-nd 


0.04 


0.23 


0.08 


0.65 



Table 6. Overview of the core-collapse supernovae which occur- 
ring in models W4 and W6. The progenitors are presented in the 
first column. When two stars are enclosed by parentheses, the 
left star explodes. Two stars enclosed by braces are the result of 
a merger or collision. The second column identifies the supernova 
type. Subsequent columns indicate how many core-collapse super- 
novae occurred in the models. The subdivisions dyn and non-dyn 
refer to models with and without dynamics. Superscript numbers 
indicate the numbers of binaries that survive the supernova. The 
rows below iVesc indicate the numbers of escapers experiencing a 
core collapse supernova after being ejected from the cluster. The 
table does not list the 3 stars in each model that collapse into 
black holes. The * indicates that one supernova occurred after a 
phase of mass transfer to a neutron star. 





Type 




W4 




W6 




Type 


dyn 


non-dyn 


dyn 


non-dyn 


Ntot: 




38 


35 


25 


26 


wr 


Ic 


1 


4 2 





2 




11 


16 


13 


12 


13 


he 


lb 


1 


2 





2 1 


(gs, ms) 


11 


1 


7 


6 


4 


(gs, gs) 


11 


2 











(gs, bh) 


11 


1 











(he, ms) 


Ic 


3 


2 


2 1 





(he, ms) 


lb 


7 1 


5 


2 


2 1 


(he, he) 


lb 


1 





1 


2 


{wd, wd} 


la 


1 


2 


1 


1 


{gs, gs} 


11 


1 











N eac 


total 


16 




6 




gs 


II 


2 




1 




{wd, wd} 


la 


11 




2 




(gs/he, ms) 


Ib/II 


2 




3* 




{ms, ms} 


II 


1 










with expected supernova rates. Table lists supernovae by 
type; Table |5] gives the evolutionary state of the explod- 
ing star and also indicates the possible presence of a binary 
companion. For comparison, we also present the expected 
number of supernovae given the same initial conditions, but 
ignoring the effects of dynamical evolution. 

The adopted Scalo (1986) initial mass function contains 



0.41 stars per thousand with m > 25 M© and 3.1 stars per 
thousand with 8 < m < 25 Mq . In a naive evolution model, 
single stars with masses ^ 25 Mq may result in a Type Ic 
supernova, while stars more massive than 8Mq produce a 
Type II supernova. For a population of 2000 single stars we 
then expect ~ 0.82 Type Ic and ~ 6.2 Type II supernovae. 
(Note that this estimate does not include binaries, and is 
therefore inapplicable to Type la and lb supernovae.) 

We compare these numbers with purely binary evo- 
lution models of Portegies Zwart & Verbunt (1996). Ac- 
cording their model AK, the Type Ic supernova rate is en- 
hanced compared to models of only single stars. In the non- 
dynamical models (W4-nd and W6-nd), Type la and lb su- 
pernovae occur without affecting the other supernova types. 
One reason for this is the fact that many Type lb super- 
nova originate from stars that were once part of a mass- 
transferring binary. The primary star in these binaries tends 
to lose its envelope in a phase of mass transfer and explode 
in a Type lb supernova. The companion star, although ini- 
tially not massive enough to experience a supernova, may 
accrete some material from the primary star. This accreted 
material may be sufficient to raise the secondary star's mass 
over the limit for experiencing a supernova. Stellar mass is, 
in a sense, recycled to produce more supernovae (see also 
Portegies Zwart & Yungelson 1999). 

The relative numbers of supernova types Ia:Ib/c:II for 
the 8 model clusters (models W4 and W6 combined) are 
0.18:0.27:0.55; in the models where stellar dynamics is ig- 
nored these ratio are 0.05:0.34:0.61. The addition of dynam- 
ical interactions to the models enhances the Type la su- 
pernova rate at the expense of the Type Ib/c and Type II 
supernova rates. The enhancement of Type-Ib/c supernova 
in model W4 is due to close (post mass transfer) binaries, 
possibly because these binaries experience more dynamical 
encounters during the shallow collapse of the cluster core 
(see the discussion in H3.H . 

3.6.2 Binary survival rate 

Only two binaries survive their first supernova. One remains 
bound because the supernova results in a black hole, which 
does not receive a velocity kick. This binary experiences a 
second phase of mass transfer (see i|3.7l for details). The 
other surviving binary is the result of a supernova in a rather 
short-period binary with a main-sequence companion. The 
resultant neutron star receives only a mild kick and the bi- 
nary remains bound. Later this binary becomes an X-ray 
source (see H3.7l for details). 

Two supernovae are triggered by collisions between car- 
bon stars, resulting in Type la events, and one Type II su- 
pernova follows a collision between two supergiants. One 
neutron star is shot into its 1.9 Mq main sequence compan- 
ion following a supernova. The resulting collision product, 
a Thorn-Zytkow star, is ejected from the star cluster before 
collapsing to a black hole. 

Only two out of 29 binaries survive the formation of a 
neutron star and one survives the formation of a black hole. 
A total of six black holes result from type Ic supernovae. 

Figure 1151 gives the time history of the supernovae in 
models W4 and W6. First black holes (bullets) are formed, 
followed by neutron stars (circles, triangles and plus signs) . 
The supernovae in the W4 models seem to occur somewhat 



12 Simon F. Portegies Zwart et al. 



earlier than in the W6 models, but the mean time at which 
a supernova occurs, 23.3 Myr and 23.5 Myr for models W4 
and W6, respectively, are not significantly different. The for- 
mation of a black hole at 46 Myr in model W6 is the result 
of a radio pulsar colliding with its carbon-star companion 
(see ET71 . 



3.7 Mass transfer and peculiar binaries 

3.7.1 First Roche-lobe contact 

The high binary fractions in our models result in frequent 
episodes of mass transfer. This is illustrated in Figure 1161 
which shows the distribution of the times of first Roche- 
lobe contact in models W4 and W6. The various symbols 
indicate the state of the donor at the moment of first contact. 
(Note that helium stars, although likely to be the remnant 
of an earlier phase of mass transfer, are still counted in this 
figure.). 

In most cases the mass is transferred to a main-sequence 
star, except in the cases indicated with x , where the mass 
transfer is onto a compact object (mostly a white dwarf). In 
the majority of these latter cases the donor is a (sub)giant, 
but there are two occasions where the donor is a helium 
star, and in one case a main-sequence star (see Figure IT7I 
One episode of mass transfer occurs between a (sub) giant 
and a black hole (model W4 at t = 829 Myr). 

The time histories of models W4 are significantly differ- 
ent from those of models W6, according to a \ 2 test on the 
cumulative distributions of Figure flUl In the W4 models, 
mass transfer tends to occur earlier than in the W6 models. 
As discussed above, the difference is a result of the greater 
dynamical activity during the early phase of model W4. In 
fact, the higher binary activity in this case is caused by the 
early phase of shallow core collapse in the W4 models, which 
does not occur in the W6 models. 



3.7.2 X-ray binaries 

Mass transfer from a hydrogen- or helium-burning star onto 
a compact object generally leads to an X-ray phase. Figure 
1171 shows the distribution of orbital periods (in days) at the 
onset of Roche-lobe overflow for binaries with white-dwarf 
(in one case a black hole) accretor. Note the clear distinction 
between main-sequence and helium-star donors at short or- 
bital periods and giant donors at larger periods. The wide 
systems with giant donors are probably most comparable to 
the supersoft sources (X-ray binaries where a white dwarf 
accretes at a super-Eddington rate for a solar- mass object, 
emitting at a rate of about 10 38 erg/s). 

A peculiar object which might be observable as an X- 
ray source results from a binary in which a helium giant 
(2.5 Mq) transfers mass to a 1.17Mq carbon-oxygen white 
dwarf. First contact occurs at 45.9 Myr at an orbital pe- 
riod of 1.7 days. When the carbon star explodes and be- 
comes a neutron star (at 46.07 Myr) it receives a kick of 
only 16 km s _1 . At that moment the companion star still 
fills its Roche lobe, and as a result the two stars coalesce 
into a single object. The merger product collapses a little 
later to form a 2.42 Mq black hole. 
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Figure 18. Orbital periods of the outer components of all hi- 
erarchical triples in all models, solid curves for the multiples in 
models W4, the dotted lines for the W6 multiples. 



3.8 Triples and higher order systems 

Long-lived triples and higher-order systems pose severe chal- 
lenges to any numerical code. We encountered a total of 41 
long-lived multiples in our runs, 21 triples [4 quadruples] 
in W4 and 13 triples [3 quadruples] in W6. Only one of 
the quadruple systems is hierarchical; the others are binary- 
binary systems. 

The first multiple systems form as early as 20 Myr (by 
a binary-binary interaction); some survive for as long as 
900 Myr. Figure [TBI presents an overview of all triples in all 
models from t = to 600 Myr. (At later times the numbers 
of triples decrease as the clusters dissolve.) The orbital pe- 
riod of the outer binary is generally rather constant. This is 
not surprising as, for the triple to be long-lived, it must be 
hierarchical and isolated, and the orbital energies and hence 
periods are adiabatic invariants. The systems showing a sig- 
nificant period derivative are driven by binary evolution or 
by temporary close encounters with other cluster members. 

The fractions of long-lived triple and quadruple systems 
in our models are only ~ 0.26% and ~ 0.12%, respectively, 
considerably smaller than the fraction of triples observed in 
the Galactic field (2.6%, Duquennoy & Mayor 1991) or in 
the Pleiades cluster (~ 2%, Mermilliod et al. 1990, 1992). 
Note also that these are the total numbers of multiple sys- 
tems created at any time during our simulations. The mean 
lifetime of such a system is 280 ± 180 Myr. 

We can study the population of short-lived multiple sys- 
tems by selecting a bound pair of stars and searching for a 
third nearby star which forms a bound subsystem. At any 
moment in time the number of such short-lived multiple sys- 
tems is about twice as high as the number of persistent 
triples. At t = 600 Myr we counted a total of 46 bound 
multiple systems in all models (W4 and W6), implying in 
a triple frequency of about 0.6% per star. For comparison, 
Kroupa (1995) found a triple fraction of 0.5% using a simi- 
lar technique. The outer orbits of these triples, however, are 
generally rather soft. 

Nine collision products in our models W4 became mem- 
bers of persistent triples or quadruples. These collisions oc- 
cured in transient triples. In about half (5) the cases, the 
two stars in the inner binary collided while the intruding 
star remained bound to the merged star. (For definiteness, 
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Figure 15. Supernova history for models W4 (bottom) and W6 (top). Type lb, Ic and II supcrnovac arc identified with circles, triangles, 
and plus-signs, respectively. Black-hole formation is identified with •. Note the collapse of a neutron star to a black hole in one of the 
W6 models around t ~ 46Myr (filled triangle). The time histories of these two distributions are not significantly different. 
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Figure 16. Times of first mass transfer for binaries in models W4 and W6. The circle size indicates the evolutionary state of the donor, 
from main-sequence (small) via Hertzsprung gap, and sub-giant, with supergiant indicated by the largest circles. The X symbols indicate 
mass transfer onto a white dwarf (or a black hole, for the second X from the right in model W4). 



we take the "inner" orbit in a quadruple system to be the 
larger of the two inner orbits.) In a few cases (3), the merger 
product and its companion became bound to another third 
star, forming a stable triple. In two of these cases, the en- 
counter resulted in second collision with the merger product. 
Both these triple collisions resulted in blue stragglers which 
had, at some time, a mass exceeding twice the turnoff mass 
of the cluster. 

Of the total of 117 collisions in the W4 models, nine oc- 
cured in transient triple systems. With a triple frequency of 
only ~ 0.26%, one naively expects less than one collision to 
occur in triples, making the frequency of collisions in triples 
unusually high. In part, this higher fraction of collisions in 
multiples is due to the dynamics — such systems generally 
form close to the cluster center, where the stellar density is 
highest and most collisions are expected. In addition, col- 
lision products tend to be found near the center, because 
they are more massive than average. The higher combined 
mass of three stars and the relatively large cross section of 
the outer orbit also increases the collision probability. Sim- 
ilar circumstances led to the high collision rates observed 
in the R136 simulations of Portegies Zwart et al. (1999, see 
also Portegies Zwart & McMillan 2002). Thus, triples are 
favorable for collisions, and collision products favor triple 
systems. 

We can characterize the triple systems formed in our 
simulations as follows: All inner orbits have periods less then 



100 year with an average of about 5.1 years. All outer orbital 
periods exceed 100 year and have an average of about 4600 
years. The eccentricities of the inner and outer orbits are 
quite distinct. The inner orbits have generally rather small 
eccentricies (e = 0.44 ±0.3) and its distribution is consistent 
at a 98.8% level with the inner eccentricity distribution of 
observed triples (Tokovinin 1997). The distribution of outer 
eccentricities (e = 0.65 ± 0.14) is consistent with thermal. 



4 COMPARISON WITH PREVIOUS WORK 
4.1 McMillan & Hut 1994 

The simulations reported by McMillan and Hut (1994) used 
up to 2048 stars, with up to 20% (rather soft: 1 — 20 or 
1 — 100/cT) primordial binaries, and included the Galac- 
tic tidal field. However, they excluded stellar evolution and 
hence any stellar mass loss. In the absence of a physical 
time scale associated with stellar evolution, they presented 
their results in units of the initial relaxation time. Our W6 
models have half lives of about 6 initial relaxation times, 
much shorter than the ~ 12-30 initial relaxation times for 
the most comparable McMillan & Hut models. Stellar mass 
loss is the main reason for the more rapid dissolution of 
our models; in addition, the McMillan & Hut runs all began 
with the cluster well inside its Jacobi surface, so the clusters 
had to expand significantly before significant tidal mass loss 
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Figure 17. Orbital period (in days) at the onset of mass transfer onto a compact object, for models W4 and W6. The circle size indicate 
the evolutionary state of the donor, from main sequence (smallest circles) to supergiant (large circles). The X indicates that the donor 
is a helium star; the filled triangle indicates the period of the simgle black hole X-ray binary. 



occurred. As discussed previously, all the McMillan & Hut 
models experienced core collapse, which is not seen in our 
simulations, where core collapse is arrested by stellar mass 
loss. 

McMillan & Hut found that the binary fraction in their 
simulations first fell due to as binaries were destroyed by 
interactions with other binaries in the cluster core, then 
rose again at late times as the cluster evaporated in the 
Galactic tidal field. In contrast, we find that the (initially 
~ 50%) binary fraction in our models remains roughly con- 
stant throughout the calculation, then increases significantly 
to ~ 70% when only ~ 10% of the cluster mass remains. 
These differences are most likely attributable to the lack of 
significant core collapse and the stronger tidal fields in our 
runs. 

One of the most interesting conclusions made by McMil- 
lan & Hut is that the spatial distribution of hard bina- 
ries is different than the distribution of soft binaries (see 
their Figure 9). In Figure |7| we showed that hard binaries 
E > 1000/cT) are slightly more centrally concentrated than 
average. 

4.2 Kroupa 1995 

The main difference between our calculations and those of 
Kroupa (1995) is his neglect of binary evolution and his 
inclusion of a prescription for pre-main-sequence binaries. 
The eccentricity distribution of his 'primordial' binaries was 
therefore somewhat different from ours, but this affects only 
binaries with the highest eccentricities and the shortest or- 
bital periods. On the other hand, his neglect of binary evolu- 
tion underestimated the fraction of tidally circularized bina- 
ries. In his calculations the boundary between hard and soft 
binaries was at an orbital period of roughly 2.7 x 10 4 years, 
similar to that in our models W4. We have already compared 
the characteristics of the populations of triple systems in our 
calculations (sec H3.8I . 

4.3 De La Fuente Marcos 1997 

De la Fuente Marcos (1997) performed studies of the evolu- 
tion of small (TV ^ 750), tidally limited open clusters hav- 
ing a variety of initial mass functions, with and without the 
inclusion of stellar (but not binary) evolution. All models 
started with a substantial fraction (1/3) of primordial bina- 
ries having mass ratios of 0.5 and energies in the ~ l-10fcT 
range. He found that the dissolution time scale of his mod- 
els depended quite sensitively on the choice of IMF, and 
that the binary population shortly before dissolution could 



show characteristic features allowing remnants of rich and 
poor clusters to be distinguished observationally. It is un- 
clear how these results extend to larger systems. 

The binary fraction in the de la Fuente Marcos simula- 
tions ranged from 33% initially to about 51% ± 0.19% near 
the disruption of the cluster (averaged over all 20 of his 
simulations). By the time the clusters dissolved, the binary 
fraction in the core had dropped to 15% ± 0.13%. We find 
a similar increase in the total binary fraction, but clearly 
the core binary fraction in our models is much higher (see 
Figure . The small binary fractions in the core near the 
end of his simulations are the direct result of his choice of 
initial binary binding energies. 

4.4 Papers I, II 

The comparison between the relative collision frequencies in 
Tab.^Jand those of model S in Paper I are quite striking. In 
that paper only encounters between single stars were stud- 
ied and primordial binaries were neglected. Still, in Paper 
I 77% of all collisions occurred between two main-sequence 
stars, and 14% between a main-sequence star and a giant; 
the remaining 9% involved giants and remnants. 

The (ms, ms) merger rates in our models are compara- 
ble to the merger rates in paper II, but (ms, wd) were much 
greater in model S of paper II because it is older (lOGyr). 
The (ms, gs) merger rate is lower because there are fraction- 
ally fewer main-sequence stars in paper II. 

4.5 Hurley et al. 2001 

Hurley et al. (2001) have modeled the open cluster M67, 
using a simulation code similar in many ways to our own. 
Both codes use fitting formulae for single star evolution, 
and include a series of recipes for dealing with binary star 
evolution. The main difference between Hurley et al. (2001) 
and the current paper lies in the parameter choices for their 
simulations: they evolved 15,000 stars to an age of 2500 Myr, 
before adding dynamics, in the form of a direct TV-body 
calculation. They then evolved the combined model to an 
age of 4300 Myr, appropriate for the old open cluster M67. 
They also present the results from smaller TV-body runs, 
starting at earlier times. 

Hurley et al. (2001) focused mostly on the formation of 
blue stragglers, and reported that their simulations indeed 
produced roughly the right number of blue stragglers, in 
agreement with observations of M67. They also performed 
a non-dynamical binary evolution population synthesis, and 
found that the number of blue stragglers in that case fell 
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short of that required by observations. They concluded that 
dynamical encounters have been crucial in the evolution of 
M67. While it is difficult to make a quantitative comparison 
between their simulations and ours, given the rather differ- 
ent types of initial conditions, our main results concerning 
the numbers and types of collisions, and the formation rate 
of triples are in broad agreement. 



5 SUMMARY AND CONCLUSIONS 

We have performed detailed Af-body calculations of 
intermediate-mass open clusters near the Sun. The initial 
conditions were selected to mimic star clusters such as 
Pleiades, Praesepe and Hyades. Our calculations included 
the effects of dynamical encounters between stars and higher 
order systems, the tidal field of the parent Galaxy and the 
evolution of single stars and binary stars. 

5.1 Cluster lifetime and structure 

Our model star clusters dissolved in the tidal field of the 
Galaxy within about f billion years. The rate of mass loss re- 
mained roughly constant, at ~ 0.8-1.4 Mq per million years, 
corresponding to a cluster half-life of ~ 600-1000 Myr, de- 
pending on distance from the Galactic center. 

The density profiles of the model clusters changed dra- 
matically during the cluster lifetime. Core collapse was pre- 
vented by stellar mass loss and binary heating. Our less con- 
centrated (W4) models became more compact during the 
first few million years, whereas the more concentrated (W6) 
models expanded from the beginning. The W4 models still 
dissolved more quickly in the Galactic tidal field, however, 
mainly due to their closer proximity of the Galactic center. 

Mass segregation is a very efficient process. After only 
a small fraction of an initial relaxation time, single white 
dwarfs, giants, collision products, hard binaries, mass trans- 
ferring binaries and binaries with one or two massive com- 
ponents (compared to the mean mass in the cluster) were 
all noticeably more centrally concentrated than the average 
star, as measured by the half mass radii of the various com- 
ponents. 

5.2 Binary stars 

The numbers of binaries decreased with time at roughly the 
same rate as the numbers of single stars, with the result that 
the binary fractions in our models remained more or less con- 
stant over the studied range of cluster ages. The binary frac- 
tion generally decreased during the first few million years by 
several percent, due to supernovae and dynamical encoun- 
ters, both of which tended to disrupt binaries. For the re- 
mainder of the evolution, this fraction slowly increased. The 
maximum binary fraction of ~ 65% is reached near disrup- 
tion. The fraction of binaries in the core remained between 
50% and 60% higher than the fraction near the half-mass 
radius over the entire lifetime of the cluster. 

The widest binaries were dissociated soon after the start 
of the simulations. A considerable fraction (4.5 ± 0.6%) of 
highly eccentric binaries with rather short orbital periods 
circularized and some experienced mass transfer early in the 
evolution of the cluster. These short-period binaries were 



generally more centrally concentrated than wider binaries 
or single stars, because they are on average more massive 
than single stars and do not interact. 

The overall distributions of binary parameters, however, 
hardly change with time. Apart from the loss of binaries with 
the largest orbital periods and those with short orbital peri- 
ods and high eccentricities, most binaries were relatively un- 
affected by either stellar evolution or stellar dynamics. The 
observed binary populations in open clusters may therefore 
be a reasonable representation of the initial primordial pop- 
ulation. 

5.3 Escapers 

We define three families of escaping stars. (1) Neutron stars 
were ejected from the cluster with very high velocities. 
About 0.4% of all stars were ejected as neutron stars due 
to supernova explosions. (2) Massive stars tended to have 
rather high escape velocities, because they were often found 
in close binaries which were disrupted by the explosion of the 
primary star, possibly after a phase of mass transfer. The 
average escape velocity of all single main sequence stars in 
model W4 was v CBC — 1.2 ± 0.51 km s _1 , whereas stars with 
m > 4M had v csc = 17 A ± 2.3 km s _1 . For model W6 
these numbers are only slightly smaller. (3) The mean es- 
caper velocity of low-mass stars (v CBC = 1.0 ± 0.28 km s _1 
for stars with m < 0.5 Mq) was comparable to the clus- 
ter velocity distribution; the mean velocity of escaping bi- 
naries was (iW = 0.97 ± 0.29 km s" 1 for model W4 and 
v CBC = 0.81 ±0.29 km s _1 for model W6), somewhat smaller 
than that for single stars. 

5.4 Stellar collisions 

Not surprisingly, collisions tended to occur near the clus- 
ter center, but a considerable fraction of mergers occured 
farther out, near the tidal radius. More than 80% of all col- 
lisions occurred within about two core radii. Most (~ 80%) 
collisions occurred between two main sequence stars. About 
70% of the collisions were the result of an unstable phase 
of mass transfer. The rest occurred in single-star-binary or 
binary-binary encounters. Multiple collisions were rare: only 
two out of 241 in our model calculations. 

Collision participants were somewhat more massive 
than expected based on the simple cross section arguments 
presented in Paper III. Most (85%) collisions occurred be- 
tween main-sequence stars; most of the remainder involved 
a main-sequence star and a giant (12%). Only 3% of the col- 
lisions occurred between two remnants. The dynamical ac- 
tivity in the cluster core, however, caused quite a few white 
dwarfs to become binary members, which would later merge 
due to the emission of gravitational waves. This leads to a 
rather high merger rate among white dwarfs (see also Shara 
& Hurley 2002). 

5.5 Supernovae and stellar remnants 

The overall supernova rate in our simulations was consis- 
tent with rates observed in the Galaxy. However, the ratios 
of Type II to Type la, lb, and Ic supernova are quite differ- 
ent. Most striking is the order of magnitude enhancement 
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of type la supernovae in the W4 models compared to the 
non-dynamical models and the population synthesis of non- 
interacting binaries. In these same models, Type lb super- 
nova are enhanced by about a factor of two compared to a 
population of binaries without dynamical encounters. 

Because of the assumed neutron-star kick distribution, 
pulsars were generally not retained in our cluster models. 
Only one out of 61 neutron stars was retained to form an 
X-ray binary later in its lifetime. The observed velocity dis- 
tribution of ejected pulsars was somewhat smaller than the 
input kick velocity distribution. X-ray binaries in which a 
white dwarf accretes from a companion star were about an 
order of magnitude more common. X-ray binaries contain- 
ing black holes were rare simply because such objects form 
from the highest-mass stars, which are themselves rare. 

Among the more unusual collisions, in our simulations, 
we observed one between a neutron star and a helium giant, 
two between CO white dwarfs, resulting in type la super- 
novae and the formation of an X-ray binary with a black 
hole as accretor. Several white dwarfs experienced a phase 
of mass transfer from a companion star. The donor was most 
often a (sub)giant. In total, 13 X-ray binaries were formed, 
most with a white dwarf as accretor. 



5.6 Multiple systems 

The number of persistent triples formed in our calculations 
was about an order of magnitude smaller than the observed 
fraction of triples in open clusters and the Galactic disk, 
although the frequency of triples and higher-order systems 
present at any moment in our simulations was comparable 
to the numbers actually observed. The majority of these 
were transient, however, and their orbital parameters dif- 
fered considerably from observations. The orbital param- 
eters of the outer orbits of the formed triples formed in 
our models was not representative of triples observed in the 
Galaxy, which tend to have shorter periods and lower eccen- 
tricities. 
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ABSTRACT 

We have modeled in detail the evolution of rich open star clusters such as the Pleiades, 
Praesepe and Hyades, using simulations that include stellar dynamics as well as the 
effects of stellar evolution. The dynamics is modeled via direct TV-body integration, 
while the evolution of single stars and binaries is followed through the use of fitting 
formulae and recipes. The feedback of stellar and binary evolution on the dynamical 
evolution of the stellar system is taken into account self-consistently. 

Our model clusters dissolve in the tidal field of the Galaxy in a time span on the 
order of a billion years. The rate of mass loss is rather constant, 1 ~ 2 M Q per million 
years. The binary fraction at first is nearly constant in time, then increases slowly near 
the end of a cluster's lifetime. For clusters which are more than about 10 8 years old 
the fractions of stars in the form of binaries, giants and collision products in the inner 
few core radii are considerably higher than in the outer regions, beyond the cluster's 
half mass radius. When stars with masses ^ 2 M Q escape from the cluster, they tend 
to do so with velocities higher than average. 

The stellar merger rate in our models is roughly one per 30 million years. Most 
mergers are the result of unstable mass transfer in close binaries (~ 70%), but a 
significant minority are caused by direct encounters between single and binary stars. 
While most collisions occur within the cluster core, even beyond the half mass radius 
collisions occasionally take place. We notice a significant birthrate of X-ray binaries, 
most containing a white dwarf as the donor. We also find some X-ray binaries with a 
neutron-star donor, but they are relatively rare. The persistent triple and higher order 
systems formed in our models by dynamical encounters between binaries and single 
stars are not representative for the multiple systems observed in the Galactic disk 
or in open clusters. We conclude that the majority of multiples in the disk probably 
formed when the stars were born, rather than through later dynamical interactions. 

Key words: Methods: N-body simulations - Binaries: general - Binaries: close - 
Open cluster and associations: general - Open cluster and associations: NGC2516, 
NGC2287, Praesepe, Hyades, NGC 2660, NGC 3680 



1 INTRODUCTION 

Open clusters are useful laboratories for studying the inter- 
play between single star evolution, binary star evolution and 
stellar dynamics. Unlike their bigger (and older) siblings, 
the globular clusters, they contain a manageable number of 
stars, and the evolution of the majority of the stars is not 
expected to be strongly affected by the dynamics. Still, a 
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significant number of collisions and subsequent stellar merg- 
ers can take place, as well as dynamically induced exchange 
reactions between single stars and binaries. For all these 
reasons, the usual zoo of objects created in binary stellar 
evolution is significantly enlarged by the presence of even 
more exotic specimens that could not have formed in vitro 
through isolated evolution, but only in vivo through the dy- 
namical interplay of initially unrelated (single or multiple) 
stars. 

The simulations reported here have been run on a 
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Table 1. Observed and derived parameters for several open star clusters with which our simulations may be compared. References to 
the literature (second column) are: (a) Pinfield et al. (1998); Raboud & Mermilliod, (1998); Bouvier et al (1998); (b) Abt & Levy (1972); 
Dachs, J & Kabus (1989); Hawley et al. (1999). (c) Harris et al. (1993); Ianna et al (1987); Cox (1954). (d) Andrievsky (1998); Jones & 
Stauffer (1991); Mermilliod & Mayor (1999); Mermilliod et al. (1990); Hodgkin et al. (1999). (e) Perryman et al. (1998 and references 
therein) Reid & Hawley (1999); (f) Frandsen et al. (1989); Hartwick & Hesser (1971); Sandrelli et al. (1999). (g) Hawley et al. 1999 ; 
Nordstrom et al. (1997); Nordstrom et al. (1996), Subsequent columns give (3) the distance to the cluster (in pc), (4) the cluster age (in 
Myr), (5) the half mass relaxation time, (6) the crossing time, (7) the total mass (in Mq), (8) estimate for the half mass radius (in pc), 
and (9) the core radius. In cases where the relaxation time is not given in the literature, we calculate it (see Paper IV); these entries are 
indicated by *. The final two columns contain information about cluster stellar content. The column labeled / s: b:t indicates the number 
of single stars, binaries and triples (separated by colons). For clusters where the numbers are given directly by observations, the table 
gives the observed numbers of each system. If the binary fraction is derived by other methods, we give the relative fractions normalized 
to the number of single stars. The last column (-/V bss:Be:gs ) gives the number of observe blue stragglers, Be stars and giants, separated 
by a column. 



name ref. d t t rix t hm M r tide r hm r CO re / s: b:t N bss:B e:g S 

[pc] -[Myr]- [M ] - [pc] - 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) 



Pleiades 


a 


135 


115 


150 


8 ~1500 


16 


2-4 


1.4 


137:60:2 


3:3:? 


NGC2516 


b 


373 


110 


220 


1000 




2.9 




16:6:? 


6:3:4 


NGC2287 


c 


675 


160-200 






6.3 






1:0.6:? 


3:1:8 


Praesepe 


d 


174 


400-900 


220 


1160 




12 


2.8 


1:0.3:0.03 


?:5:? 


Hyades 


e 


46 


625 


320 


15 1027 


10.3 


3.7 1 


2.6 


1:0.4:0 


1:0:4 


NGC 2660 


f 


2884 


900-1200 


260 


400 


? 


4 


~ 1 


0:0.3:? 


18:4:39 


NGC 3680 


g 


735 


1450 


28 


100 


4.3 


1.2 


~ 1 


44:25:0 


4:4:17 



Table 2. Initial conditions and parameters for the selected models. The first column gives the model name, followed by the cluster mass 
(in Mq), the King (1966) parameter Wo, the distance to the Galactic center (in kpc), the initial relaxation time and half mass crossing 
time (both in Myr). The remaining columns give the location of the cluster's Jacobi surface along the X— (towards the Galactic center), 
Y— and Z— (towards the Galactic pole) axes, and the initial virial radius, half mass radius and core radius (all in parsec). 



name 


M 


W 




*rlx *hm 




'"Jacob 








''core 




[M ] 




[kpc] 


[Myr] 




[PC] 






[PC] 




W4 


1708 


4 


6.8 


135 4.07 


14.5 


9.7 


7.2 


2.5 


2.14 


0.83 


W6 


1603 


6 


10.4 


140 4.15 


21.6 


14.4 


10.8 


2.5 


2.00 


0.59 



GRAPE-4 (Makino et al. 1997) system, a special-purpose 
computer designed to speed up stellar-dynamical calcula- 
tions. While the models in this paper contain only ~ 3,000 
stars, appropriate for open clusters, we have started to use 
the next-generation special-purpose computer, the GRAPE- 
6, to extend our simulations to include one or two orders of 
magnitude more stars. This will enable us to model glob- 
ular clusters on a star-by-star basis, including those with 
very dense cores where most stars are influenced by various 
'traffic accidents' in the form of encounters and mergers. 
The simulations reported here thus play a double role: as- 
trophysically, they provide insight in the evolution of open 
clusters, and computationally, they help pave the way for 
the modeling of the more complex globular clusters. 

In this paper we describe a series of self-consistent sim- 
ulations of the evolution of open star clusters. All important 
effects are included to some degree of realism: stellar evolu- 
tion, binary evolution, the internal dynamical evolution and 
the effects of the tidal field of the Galaxy. This is the fifth 
paper in a series in which we have gradually increased the 
'ecological' complexity of stellar interactions to a realistic 
level. In this paper we analyze the results of the same eight 
calculations performed for our previous Paper IV (Porte- 
gies Zwart et al. 1999). There we concentrated on global 



cluster properties and compared them with observed open 
clusters, such as the Pleiades, Praesepe and Hyades. Specifi- 
cally, we took a photometric standpoint and studied changes 
in the Hertzsprung-Russell diagram and cluster morphology 
as functions of time and initial conditions. Here we concen- 
trate on the binary population and on higher-order multiple 
systems, and we study the dynamical and observational ef- 
fects of these binaries and multiples on the evolution of the 
stellar system as a whole. 

Detailed descriptions of the numerical methods used 
and global assumptions made in our calculations are given 
in Paper IV. Since we analyze the same data from a different 
perspective we do not repeat this information here. Instead, 
we quickly review the initial conditions of our models, and 
restate the methods used; for more details, see Paper IV. 



2 METHODS 

In this section we briefly discuss the initial conditions of 
our models and the numerical techniques used in our model 
calculations. For details, refer to Appendix B and C of Paper 
IV (Portegies Zwart et al. 1999). 
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Table 3. Initial conditions for the stellar and binary population. 
The first column gives the parameter, the second column gives 
the functional dependence, followed by the lower and upper limits 
adopted. (RLOF indicates that the initial binary may be semi- 
detached.) 



parameter 


function 


lower 


upper 








limits 


mass function 


Scalo 


0.1 M Q 


100 M 


secondary mass 


P(m) = constant 


0.1 Mq 


M prim /M 


orbital separation 


P(a) = 1/a 


RLOF 


5000 AU 


eccentricity 


P(e) = 2e 





1 



solar neighborhood, we find that a model with Wo = 6 has 
to be placed slightly farther (10.4 kpc) from the Galactic 
center than the Sun, while a model with Wo = 4 has to 
reside somewhat closer (6.8 kpc), in order to obey the above 
relationships. 

For a total cluster mass of 1600 Mq, the Lagrange points 
of our two standard clusters lie at distances 14.5 pc (Wo = 
4) and 21.6 pc (Wo = 6), from the cluster center. Stars 
are removed from a simulation when their distance to the 
cluster's density center exceeds twice the distance to the first 
Lagrangian point. Table 2 reviews the selected parameters 
and initial conditions. 



2.1 Initial conditions 

We are interested in moderately rich (~ 1000 stars) open 
clusters of intermediate age ( £ lGyr). Starting from the 
currently observed mass and dynamical properties of such 
clusters we have reconstructed two plausible sets of initial 
conditions, which are presented in Table 2. Each calculation 
is performed four times with different random seeds in order 
to improve statistics. The notation in this paper is identi- 
cal to that used in Paper IV (see its Appendix A for an 
overview) . 

For our simulations we assume a Scalo (1986) initial 
mass function, with minimum and maximum masses of 
0.1 Mq and 100 Mq, respectively, and mean mass (m) ~ 
0.6Mq at birth. Consistent with the above mass estimates, 
our simulations are performed with 1024 single stars and 
1024 binaries, for a total of 3096 (3k) stars. 

Stars and binaries within our model are initialized as 
follows. A total of 2k single stars are selected from the initial 
mass function and placed in an equilibrium configuration in 
the selected density distribution (see below). We then ran- 
domly select lk stars and add a second companion star to 
them. The masses of the companions are selected randomly 
from the IMF between 0.1 Mq and the primary mass (see 
Tab. 3). For the adopted binary fraction, the restricted sec- 
ondary mass range translates into an overall mean mass of 
0.53Mq. Then the other binary parameters are determined. 
Binary eccentricities are selected from a thermal distribution 
between and 1. Orbital separations a are selected with uni- 
form probability in log a between Roche-lobe contact and an 
upper limit of 5,000 A.U. (= 10 6 R Q , or 0.02 pc). When a 
binary appears to be in contact at pericenter, that particu- 
lar orbit choice is rejected and new orbital parameters are 
selected. Table 3 gives an overview of the various distribu- 
tion functions from which stars and binaries are initialized. 
(Figure 4 shows the initial distributions of orbital period 
and eccentricity.) 

We select initial density profiles from the density dis- 
tributions described by Heggie & Ramamani (1995) with 
Wo = 4 and Wo = 6, and refer to these models as W4 
and W6, respectively, throughout this paper. We have per- 
formed four independent simulations for each set of initial 
conditions, labeled I to IV. 

All our cluster models start with the same virial radius 
of Ro = 2.5 pc. This implies a conveniently constant scaling 
between the cluster dynamical time scale ~ (GMq/Rq^ 
(= 1.5 Myr for M = 1600M Q ) and the time scale for stellar 
evolution. Each cluster is assumed to precisely fill its Jacobi 
surface at birth. Given the observed Oort constants in the 



2.2 The iV-body integrator 

The JV-body integration in our simulations is carried out 
using the kira integrator, operating within the Starlab 
software environment (McMillan & Hut 1996; Portegies 
Zwart et al. 1999). Kira uses a fourth-order Hermite scheme 
(Makino & Aarseth 1992), incorporates block time steps 
(McMillan 1986a; 1986b; Makino 1991), includes special 
treatments of close two-body and multiple encounters of ar- 
bitrary complexity, and a robust treatment of stellar and 
binary evolution and stellar collisions (see below). As men- 
tioned above, the special-purpose GRAPE-4 (Makino et al. 
1997) system is used to accelerate the computation of grav- 
itational forces between stars. 

2.3 Stellar evolution 

The evolution of the stars is adapted from the prescription 
by Portegies Zwart & Verbunt (1996, §2.1). However, some 
changes are made to the mass loss in the main-sequence 
stage for massive stars and for mass-transfer remnants. We 
follow the details of model B from Portegies Zwart & Yun- 
gelson (1998). For our treatment of stellar mass loss, see 
Portegies Zwart et al. (1998). 

Neutron stars receive a high-velocity kick at the mo- 
ment of their formation. The reasons for the occurrence of 
these high kick velocities have been discussed in detail by 
Portegies Zwart & van den Heuvel (1999) and the implica- 
tions for binary evolution are discussed by Portegies Zwart 
(2000). The distribution from which the kick velocity should 
be taken is less certain. We adopt the velocity distribution 
proposed by Harmann (1997) which has a dispersion veloc- 
ity of 450 km s _1 . Each neutron star formed received a kick 
chosen from this distribution and oriented in a random di- 
rection. 



3 RESULTS 

3.1 Overall evolution of the models 

In this section we review the global properties of models 
W4 and W6, summarizing the more extensive discussion 
presented in Paper IV. Where we discuss individual mod- 
els, we focus on the same two models discussed in Paper 
IV (model W4-IV and W6-III). In other cases the data for 
several models, or even all models, are combined in order 
to improve statistics; those cases are indicated specifically 
below. 
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Figure 1. Total mass as a function of time for model W4 (dashes) 
and W6 (solid). 



Figures 1 and 2 present the overall evolution for our 
two sets of cluster initial conditions. Figure 1 shows the 
time evolution of the cluster mass. Models W4 and W6 lose 
mass at roughly constant rates of about 1.4 Mq (0.09%) and 
0.82 M (0.05%) per million years, respectively. These mass 
loss rates result in the disruption of the cluster at about 1 
Gyr for model W4 and around 2 Gyr for the more concen- 
trated model W6. The higher mass loss rate of model W4 
is mainly attributable to its closer proximity to the Galac- 
tic center and not per se to its lower concentration. These 
rates are consistent with a mass loss rate of 1.2 M per mil- 
lion years for the slightly more massive star cluster studied 
by Hurley et al. (2001; 2002). The mass-loss rates per half 
mass relaxation time derived by Portegies Zwart et al. (2001) 
for dense star clusters near the Galactic center are about a 
factor of four higher than the corresponding rates for the 
clusters studied here. 

Figure 2 shows the evolution of the Lagrangian radii 
containing 5%, 25%, 50% and 75% of the cluster mass, 
for models W4 and W6. Both clusters expand by about 
a factor of three during the first half-mass relaxation time 
(~ 140 Myr). It is not clear whether the clusters experience 
any significant core collapse during this early phase. In Pa- 
per IV we concluded that the clusters do not experience core 
collapse due to the effects of binary activity, which tends to 
counteract core contraction (see McMillan et al. 1990, 1991). 
At later times (t ^ 650 Myr for model W4 and somewhat 
later for model W6), the Lagrangian radii decrease again as 
the cluster starts to dissolves in the tidal field of the Galaxy. 



3.2 Global binary properties 

Both models start with 50% binaries, so two-thirds of the 
stars are initially members of binary systems. The adopted 
upper limit (lO 6 !?,©) on the initial semi-major axis means 
that most initial binaries are relatively wide. The fraction of 
hard (\E\ ^ lkT) binaries was 0.46 ±0.01 for all models, ir- 
respective of the initial density profile. (The thermodynamic 
unit kT is defined in Paper IV). 

Figure 3 shows the evolution of the binary fraction. 
Shortly after the start of the runs, the binary fractions drop 
to about 42% and 48% for models W4 and W6, respectively, 
and remain roughly constant thereafter. This initial decrease 




500 750 



Figure 2. Evolution of the 5%, 25%, 50% and 75% Lagrangian 
radii for model W6 (solid lines for the 5% and 50%, dashed lines 
for the 25% and 75% Lagrangian radii). Lagrangian radii for 
model W4 are presented as dotted lines. 



+ 0A6 

y.44 

0.42 



250 



500 
t [Myr] 



750 



1000 



Figure 3. Binary fractions as a function of time for models W4 
(dotted line) and W6 (solid line). 



in binary frequency is consistent with the disruption of all 
soft (|-B| ;$ lkT) binaries. Supernova explosions dissociated 
1.8% of the binaries in the W4 models, and about half that 
number in the W6 models. This difference is mainly the 
result of random fluctuations. Fluctuations in the binary 
fraction after ~ 100 Myr are mainly the result of escapers 
and mergers. The number of mergers resulting from unstable 
mass transfer or direct collisions is 12.1 ± 0.5% for models 
W4 and W6. The escape rate of binaries closely follows the 
escape rate of single stars; both are relatively constant with 
time. Model W6 loses one binary per 2Myr; the escape rate 
for model W4 is about twice as high. 



3.3 Evolution of binary parameters 

The ecological interplay between stellar (and binary) evolu- 
tion, stellar (and binary) dynamics, and the external tidal 
influence of the Galaxy transforms the initial distributions 
of stars and binaries in complex ways. We assume that ini- 
tially all binaries are distributed throughout the cluster in 
the same way as single stars. However, since binaries are on 
average more massive than single stars, they subsequently 
tend to sink to the cluster center, where superelastic encoun- 
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ters quickly modify the spatial distributions of both stars 
and binaries. 

These encounters come in many forms, from the more 
common single-star-binary and binary-binary types to the 
rare forms that involve triples and occasionally even higher- 
order multiple systems. The combined effects of such en- 
counters significantly modifies the binary distribution func- 
tions, in terms of radius and binding energy as well as stel- 
lar type. In addition, stellar evolution and isolated binary 
evolution also cause binary parameters to evolve in time, in- 
dependent of the occurrence of dynamical close encounters. 
In this subsection, we will present some of the net results of 
these complex processes. 



3.3.1 Distributions of orbital period and eccentricity 

Figure 4a shows the initial distribution of orbital period (in 
years) and eccentricity of a subset of the primordial bina- 
ries of model W6 (the initial conditions for model W4 are 
identical). Figures 4b and c show the same distributions for 
models W4 (panel b) and W6 (panel c) at 600 Myr. In the 
middle panel, for model W4, only 871 binaries remained. In 
order to facilitate visual comparison, the upper and lower 
panels also show only 871 binaries out of the larger numbers 
that could have been plotted. 

A striking feature of the three panels in Fig 4, is the lack 
of wide (porb ~ 10 4 yr) binaries at later times. These binaries 
are completely absent in model W4, while in model W6 a 
small population of wide binaries remains. This deficiency 
in wide binaries in evolved clusters is caused by ionization of 
the soft binaries by dynamical encounters (see also Figure 
3). For model W4 this process is more efficient than for 
the more concentrated models because 1) the W4 models 
experience a high density phase during their early evolution, 
and 2) the W6 models have extended outer regions with 
relatively low stellar densities because of the weaker tidal 
field. All surviving binaries with p rt> ~ 10 6 year are located 
well beyond the cluster's half mass radius, and most contain 
at least one white dwarf. In many cases, the orbital periods 
of these non-interacting binaries has increased substantially 
due to the large amount of mass lost from one (or both) of 
the component stars. 

Binaries located above and to the left of the left-most 
dotted curves experience strong tidal effects during their 
early evolution. As a result, they quickly circularize, as can 
be seen by their precipitation to the zero-eccentricity line 
at the bottom in the last two panels. The majority of these 
systems contain at least one white dwarf or helium star. Oc- 
casionally a binary may enter the 'forbidden' left-most area 
at a later time as a result of a strong dynamical encounter. 
Such a binary will either be quickly circularized or its com- 
ponents will collide and merge to form a single star (see 
§3.5). 

Figure 5 shows how the percentage of tidally circular- 
ized binaries evolves over time. Not surprisingly, the very 
hard binaries, which have high binding energy and therefore 
tight orbits, have on average undergone much stronger tidal 
interactions and therefore have had a much larger chance to 
circularize. The four percentages presented in the figure — 
hard binaries and all binaries for each of the two classes 
of model — at first change little during the evolution of the 
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Figure 4. Scatter plots for the orbital period versus the eccen- 
tricity for 871 binaries for various models and moments in time: 
a) initial distribution, b) the distribution for model W4 at an 
age of 600 Myr, c) model W6 at t — 600 Myr. Each panel is con- 
structed using 871 binaries. Contours give the 0.5, 0.25 and 0.125 
iso-density curves of the distribution in the two coordinates plot- 
ted here. The left dotted curve indicates at which orbital period 
and eccentricity a binary with a 0.1 Mq zero-age primary circu- 
larizes during the time span of our simulations. The right dotted 
curve is for a 2.26 Mq primary at zero age. 



cluster, but they increase at later times in the case of model 
W4. 

Note that in Fig 4 there are some circularized binaries 
with long orbital periods, far to the left of the right-most 
dotted line, in the middle and lower panel. These are absent 
in the upper panel, which shows the initial conditions. The 
tidal circularization of these systems occurred when one of 
the stars ascended the giant branch before the onset of mass 
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Figure 5. The relative number of binaries that are tidally circu- 
larized, as a function of time, for the models W4 (dashed lines) 
and W6 (solid lines). In each case, the lower lines present the over- 
all fraction of binaries that are circularized, while the top curves 
indicate the fraction of circularized binaries among those binaries 
that are very hard (|£ b m ~ lOOOfcT). 




ecc 

Figure 6. Cumulative distribution of the binary eccentricity. 
Solid line, initial conditions (the thermal distribution). The 
dashed and dotted line represent the eccentricity distributions 
at an age of 600 Myr of binaries with a binding energy i?bin < 
-W 3 kT and E hin < -10 4 A;T, respectively. For the latter two 
curves the data for all models W4 and W6 are combined. How- 
ever, we excluded those binaries which have been circularized by 
tidal forces. 



transfer. While in many cases this type of evolution leads to 
a later shrinking of the orbit, in some rare cases the binaries 
remain wide after a period of modest mass transfer. 

Since tidal forces drop off quickly with distance, in most 
cases they either succeed in fully circularizing a binary, or 
they do not cause much of an effect. Consequently, the eccen- 
tricity distribution of most non-circularized binaries changes 
little with time. Figure 6 illustrates this by presenting the 
cumulative distribution of all binaries with the exception of 
those with zero eccentricity. Note that this general argument 
does not hold for the hardest binaries: when we select only 
the hard (\E bin \ ~ 1000&T) or super-hard (\E bin \ ~ 10 4 &T) 
binaries, clear deviations from the initial thermal distribu- 
tion are evident. This is not surprising: a W 4 kT binary sim- 
ply does not have enough room to develop an eccentricity 
of, say, 0.9, because the semimajor axis is too small. 



3.3.2 Radial distribution of binaries 

Figure 7 shows the binary fraction as function of distance to 
the cluster center at birth and at t = 600 Myr, adopted as 
a representative snapshot of an older cluster. (Qualitatively 
similar results would have been obtained had we selected 
any other time between ~ 200 Myr and 1 Gyr.) The upper 
solid line shows the initial distribution for all binaries, re- 
gardless of their binding energy. The lower solid line presents 
the initial fraction of binaries with a binding energy of at 
least l-Ebinl > lOOOfeT. To avoid clutter, we have shown the 
statistical uncertainties in the form of error bars only for the 
bullet points; the uncertainties in the open and close trian- 
gles can be estimated by the scatter between neighboring 
points. 

The triangles in Figure 7 shows the distribution for all 
binaries at t = 600 Myr for model W6 (open) and for model 
W4 (filled). The two distributions show similar trends, ex- 
cept for a generally smaller binary fraction in model W4 (see 
also Figure 5). The binary fraction in the central portion of 
the cluster is considerably higher than it was initially (upper 
solid line), because of mass segregation. For higher r values, 
the fraction of objects that are binary remains below the 
original value, all the way to and beyond the tidal radius. 

The bullets indicate the distribution of binaries with 
|-Ebin| > 1000&T at a cluster age of 600 Myr. In other to 
improve our statistics, we have combined the data from the 
two models; the distributions are similar, with model W6 
containing ~ 8 % more hard binaries than model W4. Note 
that we find a higher proportion of hard binaries over the 
entire cluster field as compared to the initial distribution. 
This is mainly a result of dynamical evolution: hard bina- 
ries tend to get harder, and therefore some of the binaries 
contributing to the bullet points started off with a binding 
energy less that 1000&T, and thus were not counted in the 
construction of the lower solid line. 

The higher central concentration of very hard binaries, 
as compared to less hard binaries, is evident in Figure 8, 
which figure shows the cumulative distribution for single 
stars, hard |B b in| > 1000&T and super hard |B b in| > W 6 kT 
binaries. The figure also presents the cumulative radial dis- 
tribution for higher order (mostly triple) systems. The left- 
most and right-most dotted curves in Figure 8 show the 
radial distribution for model W6 at birth and at an age of 
600 Myr, respectively. Initially, both single stars and bina- 
ries follow the left-most dotted curve. At an age of 600 Myr, 
all binaries are more centrally concentrated than the single 
stars (see the figure caption). The harder the binaries, the 
more centrally concentrated their distribution has become. 
In the case of triples, the condensation toward the center is 
extreme, with a triple half-mass radius of around 0.3 pc. 

The high central concentration of the very hard |-Et,in| > 
10 kT binaries is in part due to their early history. These 
binaries are largely the product of an unstable phase of 
mass transfer. The stellar components in these binaries were 
therefore more massive than they are at t = 600 Myr, so they 
are more more strongly affected by mass segregation, causing 
them to sink to the central regions. In contrast, the mod- 
erately hard binaries with binding energies between lOOfeT 
and WkT (not shown in Figure 8) are barely more centrally 
concentrated than the single stars. The triples are strongly 
centrally concentrated. This is a the most dramatic result 
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Figure 7. Fraction of binaries N^i n /(N SS + -/V bin ) as function 
of the distance to the cluster center. The upper solid horizontal 
line shows the initial binary fraction over the entire star cluster. 
The cut off at about 10 pc is close to the location of the tidal 
radius. (Initially there were no stars beyond the tidal radius.) 
The lower solid curve presents the initial distribution for binaries 
with |-Ebin| > 1000A;T. The open triangles give the distribution of 
all binaries for model W6 at an age of 600 Myr. The solid triangles 
present the distribution of all binaries for model W4 at an age 
of 600 Myr. The bullets (with 1 a Poissonian error bars) present 
the fraction of hard |-Ebin| > lOOO&T binaries as a function of the 
distance to the cluster center at an age of 600 Myr. These data 
combine models W4 and W6. To guide the eye we have placed a 
straight (dotted) line through these bullet points. 
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Figure 8. Cumulative radial distribution of single stars and bina- 
ries in models W4 and W6 (combined). The cumulative distribu- 
tion for all single stars and binaries at birth and at t = 600 Myr 
are given by the left and right dotted curves, respectively. The 
solid curve gives the distribution for binaries with -Ebin > 10 3 A;T 
at an age of 600 Myr. The dashed curve gives the same distribu- 
tion for -Bbin > 10 6 fcT binaries. The dash-3-dotted curve gives 
the radial distribution for higher order (mostly triple) systems. 
To improve statistics we accumulated all triples for models W4 
and W6 over the time range of 550 Myr to 650 Myr. 



of dynamical interactions: triple systems were not initially 
present, and must be created dynamically through encoun- 
ters. Since the encounter rate increases sharply toward the 
cluster center, triples and higher-order systems are born, 
and often quickly destroyed again, near the very center of 
the cluster. 
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Figure 9. Cumulative distribution of the velocities of escaping 
single stars and binaries in models W4 and model W6 (combined). 
The two dotted lines show the velocity distributions of single stars 
in the models at zero age (right) and at t — 600 Myr (left). The left 
solid curve presents the distribution of the escaper velocity of all 
stars and binaries integrated over time. The right solid curve gives 
the velocity distribution of stars which escaped within the first 
half-mass relaxation time (~ 150 Myr). The two dashed curves 
give the escape speeds for stars with masses m eS c > 2 Mq (upper) 
and m eS c > 4Mq (lower). 



3.4 Escapers 

The most important mechanism by which stars escape from 
the cluster is tidal stripping. Stars are assumed to have been 
tidally stripped once they reach two Jacobi (tidal) radii from 
the cluster center; they are then removed from the simula- 
tion. Stripped stars generally leave the cluster with rela- 
tively low velocities, as illustrated in Figure 9, which shows 
the distribution of escaper speeds for models W4 and W6. 

It is convenient to draw a distinction between tidal 
evaporation and dynamical ejection of stars. Operationally, 
we distinguish between the two processes by the speed of the 
escaper; an "ejected" single star or binary leaves the cluster 
with a speed exceeding the escape velocity. That is, we in- 
fer the dynamics from the speed. Typically, this velocity is 
imparted to the escaping object during a strong interaction 
with another object in the cluster core. Because of the rela- 
tively low central densities (and hence low reaction rates) in 
our model clusters, this type of ejection is rather uncommon. 
White dwarf ejection velocities are not significantly differ- 
ent from those of main sequence stars, but giants have on 
average somewhat lower escaper velocities. This is mainly a 
consequence of the fact that giants are generally large and 
therefore cannot have very close encounters without collid- 
ing, so they are less able to pick up large escape velocities. 

Another important escape mechanism is supernova ex- 
plosions. A neutron star formed in a core-collapse supernova 
is hurled into space with a "kick" velocity that can easily ex- 
ceed several hundreds of kilometers per second, greater than 
the cluster escape speed by 1-2 order of magnitude. These 
objects are hidden in the high- velocity tail of Figure 9; how- 
ever, they are clearly visible in Figure 10 as the population 
of high-speed single objects with masses between 1 M and 
2M . 

A supernova in a binary system has several effects on 
the velocity distribution of both single and binary stars. 
The sudden mass loss in the supernova event, as well as 
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the asymmetric velocity kick, affect a binary's internal or- 
bital parameters and its center of mass velocity. When the 
binary is dissociated in the supernova, the companion to the 
exploding star is ejected with its orbital velocity, while the 
newly formed compact object receives an additional velocity 
kick. The effects of supernovae on binary systems, and on 
the velocities of the binaries and single stars which result, 
are reviewed by Portegies Zwart (2000). 

Figure 9 shows the velocity distribution of escaping 
stars and binaries from models W4 and W6. Due to the 
similarities in escape speeds between the two models, we 
have combined the data in a single plot. The escaper veloc- 
ity distribution for single stars is slightly higher than that 
of the binaries. For clarity we show only the combined dis- 
tribution of single and binary escapers in Figure 9. Escaper 
velocities are generally only slightly higher than the cluster 
velocity dispersion, and somewhat higher in the W4 models 
than in W6. The small fraction of high-velocity escapers in 
Figure 9 is due to neutron stars escaping after supernovae 
explosions. 

In addition to the overall escaper velocities for the com- 
bined models, Figure 9 also shows the distribution for stars 
which escaped within the first half-mass relaxation time 
(right solid line). This distribution has a considerably higher 
mean than the overall distribution. The figure also shows 
the velocity distributions for escaping stars more massive 
than 2 Mq (top dashed curve) and m eS c > 4 Mq (lower 
dashed). These distributions are very different from the av- 
erage escaper speed. Both models W4 and W6 experience 
high-density phases during their early evolution, within 1 
initial half-mass relaxation time. The early escapers there- 
fore have higher velocities. A similar process operates for the 
more massive escapers. The turn-off mass of a 2 M star is 
about 800 Myr, while a 4 M star lives for less than about 
140 Myr. The distribution of escaper speeds for stars which 
escaped within t r ix includes the stars with m > 4M . Still, 
the distributions are considerably different, in the sense that 
the massive stars have relatively high escape speeds. The 
reason for this discrepancy is the greater dynamical activity 
of the high-mass stars, which take part much more often in 
dynamical encounters with binaries. 

The dependence of escaper velocity on escaper mass 
is illustrated in Figure 10. Most escapers have low masses 
and low velocities. The average escaper velocity in the W4 
and W6 runs was 1.14 km s _1 and 1.00 kms, respectively. 
Higher mass stars (and binaries) have considerably higher 
velocities. This trend was first observed by Blaauw (1961, 
see also Gies & Bolton 1986). Binaries are mostly ejected by 
dynamical encounters. Supernovae are also responsible for 
numerous high velocity escapers, both from kicks and from 
binary effects, as mentioned above. The neutron stars are 
clearly visible in Figure 10 as the objects having masses of 
about 1.4M and velocities up to ~1000km s _1 . 

Two rather low-mass (m ~ 0.11 Mq and m~ 0.7 Mq) 
stars have unusually large (~ 100km s _1 and 185km s _1 , 
respectively) escape velocities. The more massive of the two 
was ejected from a tight binary at about t = 5.6 Myr when 
its companion exploded in a Type lb supernova, dissociat- 
ing the binary. The low-mass object was ejected following a 
three-body encounter at t = 223 Myr. 

Of all stars ejected from models W4 and W6, only 109 
(~0.93% out of a total of roughly 12,000) have velocities 
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Figure 10. Mass versus velocity for escaping single stars (dots) 
and binaries (plus signs) in model W6. For the binaries, the total 
mass is used. 



exceeding three times the dispersion velocity Udisp of the 
parent cluster. Of these, 46 are neutron stars. However, 57% 
of escapers with masses between 4 Mq and 11 Mq, and 75% 
of escapers more massive than 11 Mq, have velocities more 
than 3«dis P - The average escaper velocity of stars between 
4Mq and 11 Mq is 8.5 km s _1 ; stars with masses exceeding 
11 Mq have an average escaper velocity of 21 km s _1 , which 
is much higher than the cluster dispersion velocity. 

These numbers are unusually high compared to the 
number of runaways in the Galactic disk. However, as our 
criterion here we have adopted three times the cluster ve- 
locity dispersion, instead of Blaauw's (1961) criterion of 
40 km s _1 (or three times the velocity dispersion in the 
Galactic disk). The fraction of high velocity stars following 
Blaauw (1961) was between ~ 1 % for stars with m ^ 11 M , 
about 2.5% for early B stars and about 20% among O stars 
(m £ 16 Mq, see also Sone 1991). Among stars of spectral 
type A, Stone (1991) found a runaway frequency of ^ 0.3 %. 
These numbers are consistent with binary population syn- 
thesis studies of binaries in which stellar dynamical encoun- 
ters are not taken into account (Portegies Zwart 2001). 

The velocity dispersion of all stars in the cluster at 
an age of 600 Myr is w d i sp = 0.72 ± 0.23 km s" 1 . For stars 
within 1 pc (about the core radius) of the cluster center, 
Vdisp = 0.98 ± 0.34 km s _1 ; stars between 10 and 15 pc of 
the center have «di sp = 0.59 ± 0.16 km s _1 . Most striking is 
the higher velocity dispersion for stars well outside the tidal 
radius (r > 15 pc) which have «di sp = 0.75 ± 0.20 km s _1 . A 
KS test indicates that the various velocity distributions are 
significantly different. 

Drukier et al. (1998) reported a similar effect in the 
Globular cluster M15, and attributed it to tidal shocking. 
However, this process is not included in our models; we 
therefore conclude that the observed increase in the veloc- 
ity dispersion outside the tidal radius cannot simply be a 
consequence of tidal shocks. 



3.5 Collisions and Coalescence 

Stellar mergers result from unstable mass transfer in a bi- 
nary system or from direct collisions between stars. Some- 
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Table 4. Mergers occurring in models W4 and W6, compared 
with those in model calculations published previously. The first 
column identifies the two stars involved in the merger. The next 
two columns give the number of mergers in each model calcula- 
tion, followed by the fraction of each merger type (expressed as 
a percentage, combining all runs). The last two columns give the 
fractions of mergers which occurred in models S of Papers I and 
II. Mergers between two remnants in Paper II (model S) are in- 
cluded here under {wd, wd}, although a few actually involved a 
neutron star or black hole. 





W4 


W6 


total 


Paper IS 


Paper IIS 


Ntot: 


117 


124 




— [%] 




{ms, ms} 


92 


110 


84 


77 


80 


{ms, gs} 


19 


11 


12 


14 


6 


{ms, wd} 





1 





5 


11 


{ms, ns} 


1 








1 


1 


{gs, ns} 











3 





{wd, wd} 


5 





2 





2 


{wd, ns} 





2 


1 









times (unstable) mass transfer is initiated by the presence 
of a third star. 

Table 4 gives an overview of the 241 collisions occurring 
during the model calculations discussed here. A schematic 
diagram of the relative frequencies of various collision con- 
figurations in presented in Figure 11. This figure may be 
compared directly with Figure 2 of Paper I, and with Figure 
5 of Paper II. For easier comparison, we have added to Ta- 
ble 4 two columns summarizing the results of Papers I and 
II. Note, however, that the comparison with papers I and II 
is not really appropriate, as these model calculations lasted 
for 10 Gyr. According to a KS test there is no significant 
difference between the time histories of the collisions in the 
W4 and the W6 models. 

Most ( ^ 80%) mergers in models W4 and W6 occur 
between two main sequence stars (see Table 4). This was 
also the case in Papers I and II. In Paper I, however, only 
collisions between single stars were considered, and the high 
proportion of main-sequence collisions simply reflects the 
high proportion of main-sequence stars. Paper II included 
both collisions between single stars and mergers resulting 
from binary mass transfer, so the results in that case should 
compare better with models W4 and W6. The collisions re- 
ported in this paper, however, have a rather small contri- 
bution from {ms, wd} mergers compared to those in Pa- 
pers I and II. This may be due in part to the smaller ages 
( £ 1 Gyr) of our models compared with Papers I and II 
(<10Gyr). 

Although clearly limited by small number statistics, an 
apparent difference between models W4 and W6 is the num- 
ber of {wd, wd} mergers, which are much more common in 
model W4. This difference is reflected in the higher Type la 
supernova rate in model W4. Interestingly, the calculations 
in Paper II also had a high {wd, wd} merger fraction. In 
fact, if the binary population of models W4 and W6 were 
allowed to evolve after the disruption of the cluster, the rate 
of {wd, wd} mergers would be very similar to the results of 
Paper II. The cause of the high proportion of white-dwarf 
mergers in model W4 and Paper II is the larger amount of 
dynamical activity in these models, which drives more bi- 
naries into a state of mass transfer, favoring the formation 
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Figure 12. Cumulative distribution of the distance from the clus- 
ter center at which collisions occur in models W4 (dashes) and 
W6 (solid). The dotted line gives the initial distribution of stars 
in the W6 models. 

of white-dwarf binaries (see also Hurley & Shara 2002). The 
resulting Type la supernova rate is discussed in the next 
section. 

Figure 12 shows the radial distribution of collisions in 
models W4 and W6. More than 80% of all collisions occur 
within the initial half-mass radius (~ 3pc). Collisions in the 
W6 models tend to be more concentrated to the cluster core 
than in the W4 models. It is striking that some (~ 5%) 
collisions occur far ( ,> 8pc) from the cluster center. These 
collisions are induced by binary evolution, whereas mergers 
in the core are mainly caused by stellar encounters. 

The expected distribution of masses and radii of col- 
liding stars can be calculated following Portegies Zwart et 
al. (1997; see their Eq. 14). Assuming a Maxwellian velocity 
distribution with velocity dispersion v and including the ef- 
fects of gravitational focusing, the number of collisions in 
the cluster per 10 8 year may be expressed as 

-(^) ! fe) 1 fe)fe)( t ^)' ,I) 

where m is the stellar mass. This rate is averaged over all 
other masses. (Note that, based on this estimate, we expect 
no collisions during the time span of our simulations for our 
adopted cluster parameters, underscoring the importance of 
binary interactions.) 

Figure 13 shows the expected (upper panel; see paper 
III for details) and observed (lower panel) distributions of 
primary and secondary masses of stars involved in collisions 
in models W4 and W6. We excluded here the stars which 
merged as a result of unstable mass transfer. The upper 
panel is created via Eq. 1 from the initial mass function and 
the same mass-radius relation for zero-age main-sequence 
stars as was used in the model calculations. Gray shades 
indicate collision probability; darker shades correspond to 
higher values. The collisions observed in our models are pre- 
sented in the lower panel. The masses of the colliding com- 
ponents in our model calculations are on average somewhat 
higher than expected. This phenomenon was first observed 
by Portegies Zwart et al. (1999; see also Portegies Zwart & 
McMillan 2002) in simulations of young star clusters having 
high central stellar densities. Portegies Zwart et al (1999) 
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Figure 11. Schematic overview of the relative frequencies of mergers between various stellar types in models W4 and W6 (combined). 
Stellar types are denoted by ms for main sequence, gs for (sub)giants and rm for white dwarfs and neutron stars. Table 4 gives the 
numbers of occurrences. 
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Figure 13. Relative collision rates between a primary and a sec- 
ondary star as expected from Eq. 1 (upper panel), and as found 
in models W4 and W6 (lower panel). For the lower panel the 
mergers which resulted from an unstable phase of mass transfer 
are excluded; only the true collisions are taken into account. The 
shading is linear in the encounter probability in the upper panel 
and in the number of collisions in the lower panel. Darker shades 
indicate higher collision frequency. Since the probability distribu- 
tion is symmetric about the line of equal masses, only the lower 
half of each figure is displayed. 



used R136, a compact and very dense star cluster in the 
30 Doradus region, as template for their simulations. We re- 
turn to this comparison in the discussion section. 

For 240 of the 241 stellar pairs experiencing a collision, 
we were able to reconstruct the orbit before the collision oc- 
curred. Figure 14 shows the orbital parameters for the bound 
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Figure 14. Orbital periods and eccentricities of coalescing bina- 
ries in models W4 and W6. Notice the rather large population of 
circular (zero eccentricity) binaries which experience a collision. 



systems thus obtained. The single pair for which this recon- 
struction was not possible experienced a collision during a 
supernova event — the neutron star formed in the supernova 
was ejected directly into its 1.9 Mq main-sequence compan- 
ion. In 181 of the 240 collisions (~ 71%), the orbit was cir- 
cular on merger. The eccentricity distribution in the other 
bound cases is consistent with a thermal distribution. This 
result is consistent with the findings in Paper II, where for 
model S, 19% of the collisions were the result of a dynamical 
encounter. 



3.6 Supernovae 

Supernovae can dramatically affect the evolution of a star 
cluster. During a supernova, the exploding star ejects a large 
fraction of its mass at high speed. This mass quickly leaves 
the cluster, as the ejection velocity far exceeds the cluster's 
escape speed. In addition, the newborn compact object may 
also receive a kick velocity large enough to escape the clus- 
ter. Binaries are very fragile to supernovae, and each binary 
hosting a supernova is likely to be dissociated. 



3.6.1 Supernova types 

Tables 5 and 6 present overviews of the supernovae observed 
in the models discussed in this paper, and compare them 
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Table 5. Observed supernova types (Obs, from Cappellaro et al. 
1997) and supernovae occurring in a standard Scalo (1986) mass 
function (model IMF), the population synthesis calculations of 
model AK by Portegies Zwart & Verbunt (1996, model PZV-AK) 
and in models W4 and W6, and when the dynamical evolution 
of the stellar system is ignored (models W4-nd and W6-nd). The 
first column indicates the model, followed by the number of su- 
pernovae of types la, lb, Ic and II. All numbers are normalized to 
1. 



model 


la 


lb 


Ic 


II 


Obs 


0.09 


0.11 


0.18 


0.62 


IMF 






0.12 


0.88 


PZV-AK 


0.03 


0.13 


0.12 


0.72 


W4 


0.22 


0.20 


0.07 


0.50 


W4-nd 


0.06 


0.20 


0.17 


0.57 


W6 


0.10 


0.20 


0.10 


0.63 


W6-nd 


0.04 


0.23 


0.08 


0.65 



Table 6. Overview of the core-collapse supernovae which occur- 
ring in models W4 and W6. The progenitors are presented in the 
first column. When two stars are enclosed by parentheses, the 
left star explodes. Two stars enclosed by braces are the result of 
a merger or collision. The second column identifies the supernova 
type. Subsequent columns indicate how many core-collapse super- 
novae occurred in the models. The subdivisions dyn and non-dyn 
refer to models with and without dynamics. Superscript numbers 
indicate the numbers of binaries that survive the supernova. The 
rows below _/V e sc indicate the numbers of escapers experiencing a 
core collapse supernova after being ejected from the cluster. The 
table does not list the 3 stars in each model that collapse into 
black holes. The * indicates that one supernova occurred after a 
phase of mass transfer to a neutron star. 





Type 




W4 




W6 




Type 


dyn 


non-dyn 


dyn 


non-dyn 


Ntot: 




38 


35 


25 


26 


wr 


Ic 


1 


4 2 





2 


gs 


II 


16 


13 


12 


13 


he 


lb 


1 


2 





2 1 


(gs, ms) 


II 


4 


7 


6 


4 


(gs, gs) 


II 


2 











(gs, bh) 


II 


1 











(he, ms) 


Ic 


3 


2 


2 1 





(he, ms) 


lb 


7 1 


5 


2 


2 1 


(he, he) 


lb 


1 





1 


2 


{wd, wd} 


la 


1 


2 


1 


1 


{gs, gs} 


II 


1 











N esc 


total 


16 




6 




gs 


II 


2 




1 




{wd, wd} 


la 


11 




2 




(gs/he, ms) 


Ib/II 


2 




3* 




{ms, ms} 


II 


1 










with expected supernova rates. Table 5 lists supernovae by 
type; Table 6 gives the evolutionary state of the explod- 
ing star and also indicates the possible presence of a binary 
companion. For comparison, we also present the expected 
number of supernovae given the same initial conditions, but 
ignoring the effects of dynamical evolution. 

The adopted Scalo (1986) initial mass function contains 



0.41 stars per thousand with m > 25 Mq and 3.1 stars per 
thousand with 8 < m < 25 Mq. In a naive evolution model, 
single stars with masses ^ 25 M may result in a Type Ic 
supernova, while stars more massive than 8 M© produce a 
Type II supernova. For a population of 2000 single stars we 
then expect ~ 0.82 Type Ic and ~ 6.2 Type II supernovae. 
(Note that this estimate does not include binaries, and is 
therefore inapplicable to Type la and lb supernovae.) 

We compare these numbers with purely binary evo- 
lution models of Portegies Zwart & Verbunt (1996). Ac- 
cording their model AK, the Type Ic supernova rate is en- 
hanced compared to models of only single stars. In the non- 
dynamical models (W4-nd and W6-nd), Type la and lb su- 
pernovae occur without affecting the other supernova types. 
One reason for this is the fact that many Type lb super- 
nova originate from stars that were once part of a mass- 
transferring binary. The primary star in these binaries tends 
to lose its envelope in a phase of mass transfer and explode 
in a Type lb supernova. The companion star, although ini- 
tially not massive enough to experience a supernova, may 
accrete some material from the primary star. This accreted 
material may be sufficient to raise the secondary star's mass 
over the limit for experiencing a supernova. Stellar mass is, 
in a sense, recycled to produce more supernovae (see also 
Portegies Zwart & Yungelson 1999). 

The relative numbers of supernova types Ia:Ib/c:II for 
the 8 model clusters (models W4 and W6 combined) are 
0.18:0.27:0.55; in the models where stellar dynamics is ig- 
nored these ratio are 0.05:0.34:0.61. The addition of dynam- 
ical interactions to the models enhances the Type la su- 
pernova rate at the expense of the Type Ib/c and Type II 
supernova rates. The enhancement of Type-Ib/c supernova 
in model W4 is due to close (post mass transfer) binaries, 
possibly because these binaries experience more dynamical 
encounters during the shallow collapse of the cluster core 
(see the discussion in §3.1). 

3.6.2 Binary survival rate 

Only two binaries survive their first supernova. One remains 
bound because the supernova results in a black hole, which 
does not receive a velocity kick. This binary experiences a 
second phase of mass transfer (see §3.7 for details). The 
other surviving binary is the result of a supernova in a rather 
short-period binary with a main-sequence companion. The 
resultant neutron star receives only a mild kick and the bi- 
nary remains bound. Later this binary becomes an X-ray 
source (see §3.7 for details). 

Two supernovae are triggered by collisions between car- 
bon stars, resulting in Type la events, and one Type II su- 
pernova follows a collision between two supergiants. One 
neutron star is shot into its 1.9 Mq main sequence compan- 
ion following a supernova. The resulting collision product, 
a Thorn-Zytkow star, is ejected from the star cluster before 
collapsing to a black hole. 

Only two out of 29 binaries survive the formation of a 
neutron star and one survives the formation of a black hole. 
A total of six black holes result from type Ic supernovae. 

Figure 15 gives the time history of the supernovae in 
models W4 and W6. First black holes (bullets) are formed, 
followed by neutron stars (circles, triangles and plus signs). 
The supernovae in the W4 models seem to occur somewhat 
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earlier than in the W6 models, but the mean time at which 
a supernova occurs, 23.3 Myr and 23.5 Myr for models W4 
and W6, respectively, are not significantly different. The for- 
mation of a black hole at 46 Myr in model W6 is the result 
of a radio pulsar colliding with its carbon-star companion 
(see §3.7). 



3.7 Mass transfer and peculiar binaries 

3.7.1 First Roche-lobe contact 

The high binary fractions in our models result in frequent 
episodes of mass transfer. This is illustrated in Figure 16, 
which shows the distribution of the times of first Roche- 
lobe contact in models W4 and W6. The various symbols 
indicate the state of the donor at the moment of first contact. 
(Note that helium stars, although likely to be the remnant 
of an earlier phase of mass transfer, are still counted in this 
figure.). 

In most cases the mass is transferred to a main-sequence 
star, except in the cases indicated with x, where the mass 
transfer is onto a compact object (mostly a white dwarf). In 
the majority of these latter cases the donor is a (sub)giant, 
but there are two occasions where the donor is a helium 
star, and in one case a main-sequence star (see Figure 17. 
One episode of mass transfer occurs between a (sub) giant 
and a black hole (model W4 at t = 829 Myr). 

The time histories of models W4 are significantly differ- 
ent from those of models W6, according to a x 2 test on the 
cumulative distributions of Figure 16. In the W4 models, 
mass transfer tends to occur earlier than in the W6 models. 
As discussed above, the difference is a result of the greater 
dynamical activity during the early phase of model W4. In 
fact, the higher binary activity in this case is caused by the 
early phase of shallow core collapse in the W4 models, which 
does not occur in the W6 models. 



3.7.2 X-ray binaries 

Mass transfer from a hydrogen- or helium-burning star onto 
a compact object generally leads to an X-ray phase. Figure 
17 shows the distribution of orbital periods (in days) at the 
onset of Roche-lobe overflow for binaries with white-dwarf 
(in one case a black hole) accretor. Note the clear distinction 
between main-sequence and helium-star donors at short or- 
bital periods and giant donors at larger periods. The wide 
systems with giant donors are probably most comparable to 
the supersoft sources (X-ray binaries where a white dwarf 
accretes at a super-Eddington rate for a solar-mass object, 
emitting at a rate of about 10 38 erg/s). 

A peculiar object which might be observable as an X- 
ray source results from a binary in which a helium giant 
(2.5 M ) transfers mass to a 1.17 M carbon-oxygen white 
dwarf. First contact occurs at 45.9 Myr at an orbital pe- 
riod of 1.7 days. When the carbon star explodes and be- 
comes a neutron star (at 46.07 Myr) it receives a kick of 
only 16 km s _1 . At that moment the companion star still 
fills its Roche lobe, and as a result the two stars coalesce 
into a single object. The merger product collapses a little 
later to form a 2.42 Mq black hole. 
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Figure 18. Orbital periods of the outer components of all hi- 
erarchical triples in all models, solid curves for the multiples in 
models W4, the dotted lines for the W6 multiples. 



3.8 Triples and higher order systems 

Long-lived triples and higher-order systems pose severe chal- 
lenges to any numerical code. We encountered a total of 41 
long-lived multiples in our runs, 21 triples [4 quadruples] 
in W4 and 13 triples [3 quadruples] in W6. Only one of 
the quadruple systems is hierarchical; the others are binary- 
binary systems. 

The first multiple systems form as early as 20 Myr (by 
a binary-binary interaction); some survive for as long as 
900 Myr. Figure 18 presents an overview of all triples in all 
models from t = to 600 Myr. (At later times the numbers 
of triples decrease as the clusters dissolve.) The orbital pe- 
riod of the outer binary is generally rather constant. This is 
not surprising as, for the triple to be long-lived, it must be 
hierarchical and isolated, and the orbital energies and hence 
periods are adiabatic invariants. The systems showing a sig- 
nificant period derivative are driven by binary evolution or 
by temporary close encounters with other cluster members. 

The fractions of long-lived triple and quadruple systems 
in our models are only ~ 0.26% and ~ 0.12%, respectively, 
considerably smaller than the fraction of triples observed in 
the Galactic field (2.6%, Duquennoy & Mayor 1991) or in 
the Pleiades cluster (~ 2%, Mermilliod et al. 1990, 1992). 
Note also that these are the total numbers of multiple sys- 
tems created at any time during our simulations. The mean 
lifetime of such a system is 280 ± 180 Myr. 

We can study the population of short-lived multiple sys- 
tems by selecting a bound pair of stars and searching for a 
third nearby star which forms a bound subsystem. At any 
moment in time the number of such short-lived multiple sys- 
tems is about twice as high as the number of persistent 
triples. At t = 600 Myr we counted a total of 46 bound 
multiple systems in all models (W4 and W6), implying in 
a triple frequency of about 0.6% per star. For comparison, 
Kroupa (1995) found a triple fraction of 0.5% using a simi- 
lar technique. The outer orbits of these triples, however, are 
generally rather soft. 

Nine collision products in our models W4 became mem- 
bers of persistent triples or quadruples. These collisions oc- 
cured in transient triples. In about half (5) the cases, the 
two stars in the inner binary collided while the intruding 
star remained bound to the merged star. (For definiteness, 
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Figure 15. Supernova history for models W4 (bottom) and W6 (top). Type lb, Ic and II supernovae are identified with circles, triangles, 
and plus-signs, respectively. Black-hole formation is identified with •. Note the collapse of a neutron star to a black hole in one of the 
W6 models around t ~ 46Myr (filled triangle). The time histories of these two distributions are not significantly different. 
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Figure 16. Times of first mass transfer for binaries in models W4 and W6. The circle size indicates the evolutionary state of the donor, 
from main-sequence (small) via Hertzsprung gap, and sub-giant, with supergiant indicated by the largest circles. The x symbols indicate 
mass transfer onto a white dwarf (or a black hole, for the second x from the right in model W4). 



we take the "inner" orbit in a quadruple system to be the 
larger of the two inner orbits.) In a few cases (3), the merger 
product and its companion became bound to another third 
star, forming a stable triple. In two of these cases, the en- 
counter resulted in second collision with the merger product. 
Both these triple collisions resulted in blue stragglers which 
had, at some time, a mass exceeding twice the turnoff mass 
of the cluster. 

Of the total of 117 collisions in the W4 models, nine oc- 
cured in transient triple systems. With a triple frequency of 
only ~ 0.26%, one naively expects less than one collision to 
occur in triples, making the frequency of collisions in triples 
unusually high. In part, this higher fraction of collisions in 
multiples is due to the dynamics — such systems generally 
form close to the cluster center, where the stellar density is 
highest and most collisions are expected. In addition, col- 
lision products tend to be found near the center, because 
they are more massive than average. The higher combined 
mass of three stars and the relatively large cross section of 
the outer orbit also increases the collision probability. Sim- 
ilar circumstances led to the high collision rates observed 
in the R136 simulations of Portegies Zwart et al. (1999, see 
also Portegies Zwart & McMillan 2002). Thus, triples are 
favorable for collisions, and collision products favor triple 
systems. 

We can characterize the triple systems formed in our 
simulations as follows: All inner orbits have periods less then 



100 year with an average of about 5.1 years. All outer orbital 
periods exceed 100 year and have an average of about 4600 
years. The eccentricities of the inner and outer orbits are 
quite distinct. The inner orbits have generally rather small 
eccentricies (e = 0.44 ±0.3) and its distribution is consistent 
at a 98.8% level with the inner eccentricity distribution of 
observed triples (Tokovinin 1997). The distribution of outer 
eccentricities (e = 0.65 ± 0.14) is consistent with thermal. 

4 COMPARISON WITH PREVIOUS WORK 
4.1 McMillan & Hut 1994 

The simulations reported by McMillan and Hut (1994) used 
up to 2048 stars, with up to 20% (rather soft: 1 — 20 or 
1 — 100/cT) primordial binaries, and included the Galac- 
tic tidal field. However, they excluded stellar evolution and 
hence any stellar mass loss. In the absence of a physical 
time scale associated with stellar evolution, they presented 
their results in units of the initial relaxation time. Our W6 
models have half lives of about 6 initial relaxation times, 
much shorter than the ~ 12-30 initial relaxation times for 
the most comparable McMillan & Hut models. Stellar mass 
loss is the main reason for the more rapid dissolution of 
our models; in addition, the McMillan & Hut runs all began 
with the cluster well inside its Jacobi surface, so the clusters 
had to expand significantly before significant tidal mass loss 
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Figure 17. Orbital period (in days) at the onset of mass transfer onto a compact object, for models W4 and W6. The circle size indicate 
the evolutionary state of the donor, from main sequence (smallest circles) to supergiant (large circles). The x indicates that the donor 
is a helium star; the filled triangle indicates the period of the simgle black hole X-ray binary. 



occurred. As discussed previously, all the McMillan & Hut 
models experienced core collapse, which is not seen in our 
simulations, where core collapse is arrested by stellar mass 
loss. 

McMillan & Hut found that the binary fraction in their 
simulations first fell due to as binaries were destroyed by 
interactions with other binaries in the cluster core, then 
rose again at late times as the cluster evaporated in the 
Galactic tidal field. In contrast, we find that the (initially 
~ 50%) binary fraction in our models remains roughly con- 
stant throughout the calculation, then increases significantly 
to ~ 70% when only ~ 10% of the cluster mass remains. 
These differences are most likely attributable to the lack of 
significant core collapse and the stronger tidal fields in our 
runs. 

One of the most interesting conclusions made by McMil- 
lan & Hut is that the spatial distribution of hard bina- 
ries is different than the distribution of soft binaries (see 
their Figure 9). In Figure 7 we showed that hard binaries 
E > 1000/cT) are slightly more centrally concentrated than 
average. 

4.2 Kroupa 1995 

The main difference between our calculations and those of 
Kroupa (1995) is his neglect of binary evolution and his 
inclusion of a prescription for pre-main-sequence binaries. 
The eccentricity distribution of his 'primordial' binaries was 
therefore somewhat different from ours, but this affects only 
binaries with the highest eccentricities and the shortest or- 
bital periods. On the other hand, his neglect of binary evolu- 
tion underestimated the fraction of tidally circularized bina- 
ries. In his calculations the boundary between hard and soft 
binaries was at an orbital period of roughly 2.7 x 10 4 years, 
similar to that in our models W4. We have already compared 
the characteristics of the populations of triple systems in our 
calculations (see §3.8). 

4.3 De La Fuente Marcos 1997 

De la Fuente Marcos (1997) performed studies of the evolu- 
tion of small (N ^ 750), tidally limited open clusters hav- 
ing a variety of initial mass functions, with and without the 
inclusion of stellar (but not binary) evolution. All models 
started with a substantial fraction (1/3) of primordial bina- 
ries having mass ratios of 0.5 and energies in the ~ l-10feT 
range. He found that the dissolution time scale of his mod- 
els depended quite sensitively on the choice of IMF, and 
that the binary population shortly before dissolution could 



show characteristic features allowing remnants of rich and 
poor clusters to be distinguished observationally. It is un- 
clear how these results extend to larger systems. 

The binary fraction in the de la Fuente Marcos simula- 
tions ranged from 33% initially to about 51% ± 0.19% near 
the disruption of the cluster (averaged over all 20 of his 
simulations). By the time the clusters dissolved, the binary 
fraction in the core had dropped to 15% ± 0.13%. We find 
a similar increase in the total binary fraction, but clearly 
the core binary fraction in our models is much higher (see 
Figure 7). The small binary fractions in the core near the 
end of his simulations are the direct result of his choice of 
initial binary binding energies. 

4.4 Papers I, II 

The comparison between the relative collision frequencies in 
Tab. 4 and those of model S in Paper I are quite striking. In 
that paper only encounters between single stars were stud- 
ied and primordial binaries were neglected. Still, in Paper 
I 77% of all collisions occurred between two main-sequence 
stars, and 14% between a main-sequence star and a giant; 
the remaining 9% involved giants and remnants. 

The (ms, ms) merger rates in our models are compara- 
ble to the merger rates in paper II, but (ms, wd) were much 
greater in model S of paper II because it is older (lOGyr). 
The (ms, gs) merger rate is lower because there are fraction- 
ally fewer main-sequence stars in paper II. 

4.5 Hurley et al. 2001 

Hurley et al. (2001) have modeled the open cluster M67, 
using a simulation code similar in many ways to our own. 
Both codes use fitting formulae for single star evolution, 
and include a series of recipes for dealing with binary star 
evolution. The main difference between Hurley et al. (2001) 
and the current paper lies in the parameter choices for their 
simulations: they evolved 15,000 stars to an age of 2500 Myr, 
before adding dynamics, in the form of a direct TV-body 
calculation. They then evolved the combined model to an 
age of 4300 Myr, appropriate for the old open cluster M67. 
They also present the results from smaller TV-body runs, 
starting at earlier times. 

Hurley et al. (2001) focused mostly on the formation of 
blue stragglers, and reported that their simulations indeed 
produced roughly the right number of blue stragglers, in 
agreement with observations of M67. They also performed 
a non-dynamical binary evolution population synthesis, and 
found that the number of blue stragglers in that case fell 
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short of that required by observations. They concluded that 
dynamical encounters have been crucial in the evolution of 
M67. While it is difficult to make a quantitative comparison 
between their simulations and ours, given the rather differ- 
ent types of initial conditions, our main results concerning 
the numbers and types of collisions, and the formation rate 
of triples are in broad agreement. 



5 SUMMARY AND CONCLUSIONS 

We have performed detailed JV-body calculations of 
intermediate-mass open clusters near the Sun. The initial 
conditions were selected to mimic star clusters such as 
Pleiades, Praesepe and Hyades. Our calculations included 
the effects of dynamical encounters between stars and higher 
order systems, the tidal field of the parent Galaxy and the 
evolution of single stars and binary stars. 

5.1 Cluster lifetime and structure 

Our model star clusters dissolved in the tidal field of the 
Galaxy within about 1 billion years. The rate of mass loss re- 
mained roughly constant, at ~ 0.8-1.4 M© per million years, 
corresponding to a cluster half-life of ~ 600-1000 Myr, de- 
pending on distance from the Galactic center. 

The density profiles of the model clusters changed dra- 
matically during the cluster lifetime. Core collapse was pre- 
vented by stellar mass loss and binary heating. Our less con- 
centrated (W4) models became more compact during the 
first few million years, whereas the more concentrated (W6) 
models expanded from the beginning. The W4 models still 
dissolved more quickly in the Galactic tidal field, however, 
mainly due to their closer proximity of the Galactic center. 

Mass segregation is a very efficient process. After only 
a small fraction of an initial relaxation time, single white 
dwarfs, giants, collision products, hard binaries, mass trans- 
ferring binaries and binaries with one or two massive com- 
ponents (compared to the mean mass in the cluster) were 
all noticeably more centrally concentrated than the average 
star, as measured by the half mass radii of the various com- 
ponents. 

5.2 Binary stars 

The numbers of binaries decreased with time at roughly the 
same rate as the numbers of single stars, with the result that 
the binary fractions in our models remained more or less con- 
stant over the studied range of cluster ages. The binary frac- 
tion generally decreased during the first few million years by 
several percent, due to supernovae and dynamical encoun- 
ters, both of which tended to disrupt binaries. For the re- 
mainder of the evolution, this fraction slowly increased. The 
maximum binary fraction of ~ 65% is reached near disrup- 
tion. The fraction of binaries in the core remained between 
50% and 60% higher than the fraction near the half-mass 
radius over the entire lifetime of the cluster. 

The widest binaries were dissociated soon after the start 
of the simulations. A considerable fraction (4.5 ± 0.6%) of 
highly eccentric binaries with rather short orbital periods 
circularized and some experienced mass transfer early in the 
evolution of the cluster. These short-period binaries were 



generally more centrally concentrated than wider binaries 
or single stars, because they are on average more massive 
than single stars and do not interact. 

The overall distributions of binary parameters, however, 
hardly change with time. Apart from the loss of binaries with 
the largest orbital periods and those with short orbital peri- 
ods and high eccentricities, most binaries were relatively un- 
affected by either stellar evolution or stellar dynamics. The 
observed binary populations in open clusters may therefore 
be a reasonable representation of the initial primordial pop- 
ulation. 

5.3 Escapers 

We define three families of escaping stars. (1) Neutron stars 
were ejected from the cluster with very high velocities. 
About 0.4% of all stars were ejected as neutron stars due 
to supernova explosions. (2) Massive stars tended to have 
rather high escape velocities, because they were often found 
in close binaries which were disrupted by the explosion of the 
primary star, possibly after a phase of mass transfer. The 
average escape velocity of all single main sequence stars in 
model W4 was w eS c = 1.2 ± 0.51 km s -1 , whereas stars with 
m > 4M had w eS c = 17.4 ± 2.3 km s" 1 . For model W6 
these numbers are only slightly smaller. (3) The mean es- 
caper velocity of low-mass stars (v eS c = 1-0 ± 0.28 km s _1 
for stars with m < 0.5 M ) was comparable to the clus- 
ter velocity distribution; the mean velocity of escaping bi- 
naries was (vesc = 0.97 ± 0.29 km s _1 for model W4 and 
Vesc = 0.81 ±0.29 km s _1 for model W6), somewhat smaller 
than that for single stars. 

5.4 Stellar collisions 

Not surprisingly, collisions tended to occur near the clus- 
ter center, but a considerable fraction of mergers occured 
farther out, near the tidal radius. More than 80% of all col- 
lisions occurred within about two core radii. Most (~ 80%) 
collisions occurred between two main sequence stars. About 
70% of the collisions were the result of an unstable phase 
of mass transfer. The rest occurred in single-star-binary or 
binary-binary encounters. Multiple collisions were rare: only 
two out of 241 in our model calculations. 

Collision participants were somewhat more massive 
than expected based on the simple cross section arguments 
presented in Paper III. Most (85%) collisions occurred be- 
tween main-sequence stars; most of the remainder involved 
a main-sequence star and a giant (12%). Only 3% of the col- 
lisions occurred between two remnants. The dynamical ac- 
tivity in the cluster core, however, caused quite a few white 
dwarfs to become binary members, which would later merge 
due to the emission of gravitational waves. This leads to a 
rather high merger rate among white dwarfs (see also Shara 
& Hurley 2002). 

5.5 Supernovae and stellar remnants 

The overall supernova rate in our simulations was consis- 
tent with rates observed in the Galaxy. However, the ratios 
of Type II to Type la, lb, and Ic supernova are quite differ- 
ent. Most striking is the order of magnitude enhancement 
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of type la supernovae in the W4 models compared to the 
non-dynamical models and the population synthesis of non- 
interacting binaries. In these same models, Type lb super- 
nova are enhanced by about a factor of two compared to a 
population of binaries without dynamical encounters. 

Because of the assumed neutron-star kick distribution, 
pulsars were generally not retained in our cluster models. 
Only one out of 61 neutron stars was retained to form an 
X-ray binary later in its lifetime. The observed velocity dis- 
tribution of ejected pulsars was somewhat smaller than the 
input kick velocity distribution. X-ray binaries in which a 
white dwarf accretes from a companion star were about an 
order of magnitude more common. X-ray binaries contain- 
ing black holes were rare simply because such objects form 
from the highest-mass stars, which are themselves rare. 

Among the more unusual collisions, in our simulations, 
we observed one between a neutron star and a helium giant, 
two between CO white dwarfs, resulting in type la super- 
novae and the formation of an X-ray binary with a black 
hole as accretor. Several white dwarfs experienced a phase 
of mass transfer from a companion star. The donor was most 
often a (sub)giant. In total, 13 X-ray binaries were formed, 
most with a white dwarf as accretor. 

5.6 Multiple systems 

The number of persistent triples formed in our calculations 
was about an order of magnitude smaller than the observed 
fraction of triples in open clusters and the Galactic disk, 
although the frequency of triples and higher-order systems 
present at any moment in our simulations was comparable 
to the numbers actually observed. The majority of these 
were transient, however, and their orbital parameters dif- 
fered considerably from observations. The orbital param- 
eters of the outer orbits of the formed triples formed in 
our models was not representative of triples observed in the 
Galaxy, which tend to have shorter periods and lower eccen- 
tricities. 
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